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Abstract 


This  study  investigated  the  possibility  of  abort  landing  the  Space  Shuttle  at  eas  :  coast  airports 
when  launched  at  inclinations  of  5  L6  degrees  or  more.  Computer  modeling  was  used  to  characterize  both 
the  Shuttle  launch  out  of  Cape  Canaveral  and  two  methods  of  un-powered  abort  descents  jfrom  various 
points  in  the  launch  following  Solid  Rocket  Booster  (SRB)  separation. 

The  first  method  used  varying  values  of  pitch  and  roll  held  constant  to  control  the  descent.  By 
plotting  the  latitude  and  longitude  of  the  point  in  the  descent  when  the  nominal  landing  altitude  was 
achieved  against  locations  of  east  coast  airports  it  was  found  that  there  are  indeed  east  coast  abort 
opportunities  for  high  inclination  launches  out  of  Cape  Canaveral. 

The  second  method  used  a  constant  pitch  and  roll  until  the  proper  heading  angle  to  intercept  a 
desired  target  airport  was  reached  then  maneuver  to  0  degrees  roll  and  30  degrees  pitch.  These  trajectories 
were  attempted  throughout  the  launch  for  different  airports  so  that  windows  of  opportunity  could  be 
established.  It  was  shown  that  these  windows  exist  but  only  for  limited  times  ranging  fi'om  8  to  34  seconds. 
These  opportunities  may  be  expanded  if  further  studies  investigate  powered  and  optimal  control  cases. 
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1.  Background 


1. 1  Introduction 

This  chapter  will  give  a  background  of  the  existing  abort  procedures  and  launch  profile  of  the 
Space  Shuttle  Flight  System.  Future  missions  of  the  Shuttle,  will  then  be  presented  as  well  as  the  problem 
statement  and  method  of  solution  of  this  thesis. 

1 .2  Space  Shuttle  Flight  System  Launch  Profile 

The  Space  Shuttle  Flight  System  is  composed  of  the  Orbiter^  an  external  liquid  propellant  tank, 
and  two  Solid  propellant  Rocket  Boosters  (SRB’s).  The  Orbiter’s  main  engines  and  the  SRB’s  provide 
propulsion  fi*om  launch  to  approximately  120  seconds  and  an  altitude  of  50km.  The  SRB’s  are  then 
released  leaving  the  Orbiter  and  external  tank  propelled  by  the  main  engines.  The  external  tank  separates 
once  the  Orbiter  is  injected  into  the  required  ascent  trajectory  ,which  is  typically  480  seconds,  after  which 
the  Orbital  Maneuvering  System  (OMS)  completes  the  insertion  into  the  desired  orbit. 

L3  Existing  Abort  Procedures 

Failure  of  one  or  more  of  the  Shuttle’s  systems  during  ascent  may  cause  the  selection  of  an  abort 
mode  to  become  necessary.  The  type  of  abort  will  depend  on  the  system  that  failed  and  the  time  in  flight. 
The  two  basic  types  of  ascent  abort  modes  are  intact  and  contingency .  Intact  aborts  are  intended  to  safely 
return  the  Orbiter  to  a  planned  landing  zone  while  contingency  aborts  generally  result  in  a  ditch  operation 
(6:1).  The  four  existing  intact  abort  procedures,  Abort  To  Orbit,  Abort  Once  Around,  Return  To  Launch 
Site,  and  TransAtlantic  Landing  can  return  the  Orbiter  safely  to  the  ground  in  as  fast  as  35  to  90  minutes 
jfrom  launch  depending  on  the  mode  selected.  Figure  1-1 .  demonstrates  the  Shuttle  launch  profile  as  well  as 
the  abort  modes. 
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L3, 1  Abort  to  Orbit,  An  Abort  To  Orbit  or  ATO  is  used  if  the  Shuttle  still  has  sufficient 
propulsion  to  place  it  into  a  lower  orbit  then  planned  The  OMS  engines  will  then  place  it  into  a  circular 
orbit  where  it  can  deorbit  as  it  would  under  normal  circumstances  and  land  at  the  planned  site.  This  mode 
is  only  chosen  if  there  are  sufficient  working  systems  and  time  to  land  is  not  a  factor  (6:1). 

1.3.2  Abort  Once  Around.  Abort  Once  Around  (AO A)  is  used  when  vehicle  performance  is 
insufficient  to  achieve  a  feasible  orbit  or  the  OMS  doesn’t  have  enough  fuel  to  perform  the  orbit  and  deorbit 
maneuvering  necessary  for  an  ATO.  Furthermore  the  AOA  is  used  in  cases  where  it  is  necessary  to  land  as 
quick  as  possible  such  as  a  cabin  leak  or  cooling  loss.  The  maneuver  is  done  by  firing  the  OMS  once  after 
main  engine  shut  off  so  that  a  second  OMS  thrust  will  send  the  vehicle  out  of  orbit  to  land  at  an  AOA 
landing  site  such  as  White  Sands,  NM,  Edwards  AFB  ,or  Kennedy  Space  Center.  This  mode  results  in  the 
Orbiter  circling  the  Earth  once  and  landing  in  90  minutes  (6:1). 
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1,3.3  Return  to  Launch  Site,  The  RTLS  abort  mode  is  selected  when  there  is  a  system  error 
during  the  first  4  minutes  and  20  seconds  of  the  ascent  at  which  time  there  is  not  enough  main  propulsive 
system  propellant  to  return  to  the  launch  site.  This  abort  selection  is  also  preferable  if  the  Orbiter  and  crew 
need  to  be  on  the  ground  quickly.  The  RTLS  mode  can  put  the  Orbiter  on  the  ground  in  25  minutes  from 
launch  (6:1). 

The  procedure  for  a  RTLS  consists  of  three  phases;  a  powered  phase,  a  separation  phase,  and  a  glide 
phase.  The  powered  phase  begins  with  the  crew  selection  of  the  abort  at  2  minutes  20  seconds  which  is  just 
after  the  separation  of  the  solid  rocket  boosters.  The  shuttle  continues  down  range  to  dissipate  excess  main 
engine  propellant  so  that  it  has  just  enough  fuel  to  turn  around,  fly  back  towards  the  launch  site,  and  achieve 
the  proper  main  engine  cutoff  conditions  so  the  vehicle  can  glide  to  the  runway  after  external  tank 
separation.  As  it  continues  downrange  the  Orbiter/extemal  tank  makes  a  pitch  around  maneuver  so  the 
thrust  of  the  main  engines  are  used  to  null  the  downrange  velocity.  All  excess  fuel  like  the  OMS  propellant 
are  dumped  at  this  time  to  improve  the  Orbiter  weight  and  center  of  gravity.  At  20  seconds  prior  to  main 
engine  cutoff  the  separation  phase  begins  with  a  powered  pitch  down  maneuver  that  brings  the  shuttle  to  the 
required  external  tank  separation  attitude  and  pitch  rate.  A  reaction  control  system  is  then  employed  to 
ensure  that  the  Orbiter  does  not  recontact  the  external  tank  and  that  the  Orbiter  has  achieved  the  necessary 
pitch  attitude  to  begin  the  glide  phase.  Once  the  external  tank  has  been  released  the  glide  phase  begins 
which  is  handled  like  a  normal  entry. 
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L3.4  Transatlantic  Landing,  The  Transadantic  Abort  Landing  (TAL)  was  developed  for  low 
inclination  launches  to  handle  emergency  situations  that  might  occur  after  the  last  RTLS  opportunity  and 
before  an  AOA  could  be  attempted.  This  mode  selection  could  either  mean  that  there  is  insufficient  power 
to  accomplish  an  AOA  or  that  the  emergency  calls  for  the  Orbiter  to  land  immediately. 

In  a  TAL  abort  the  vehicle  would  continue  out  of  Cape  Canaveral  on  a  low  inclination  trajectory 
that  would  place  the  shuttles  downrange  ground  track  in  close  proximity  to  the  West  Coast  of  Africa  or  the 
South  West  of  Europ)e.  As  shown  in  Figure  1-2  the  Shuttle  would  continue  on  a  ballistic  trajectory  across 
the  Atlantic  to  land  at  a  predetermined  site  that  lies  near  the  nominal  ascent  ground  track.  The  site  would 
have  to  have  sufficient  runway  length ,  weather  conditions,  and  approval  from  the  State  Department  (6:1). 

The  TAL  abort  mode  sends  commands  to  steer  the  vehicle  toward  the  plane  of  the  landing  site  and 
roll  heads  up  before  main  engine  cutoff.  As  in  an  RTLS,  QMS  propellant  is  dumped  to  decrease  weight  and 
place  the  center  of  gravity  in  the  proper  position  for  control.  The  landing  is  then  handled  like  a  nominal 
entry. 

1.3,5  Future  Shuttle  Missions.  The  Space  Shuttle  is  scheduled  to  launch  50  missions  through 
August  of 2003.  As  shown  in  the  pie  chart  in  Figure  1  -3.,  82  percent  of  the  launches  will  be  flown  at  a  5 1 .6 
degree  inclination  which  makes  the  TAL  abort  mode  invalid  because  the  ground  track  does  not  intercept 
any  transatlantic  landing  sites  (6: 1 ).  This  in  turn  leaves  between  4  minutes  20  seconds  until  just  prior  to 
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mam  engine  cutoff  where  the  Orbiter  has  only  a  contingency  abort  option  rather  than  an  intact  abort  option. 
The  high  inclination  launches  do,  however,  come  significantly  close  to  the  East  Coast  of  North  America, 
namely  the  United  States  of  America. 


r 


Percent  Distribution  of  Launch  Inclinations  for  Future  Space  Shuttle 

Missions 


■39  deg 
■57  deg 
■28.45  deg 
■51 .6  deg 


Figure  1-3.  Inclination  Angles  for  Future  Shuttle  Missions 


J.  4  Problem  Statement  and  Method  of  Solution 

Since  the  majority  of  the  future  launches  out  of  Cape  Canaveral  will  take  the  shuttle  close  to  the 
coast  of  the  United  States  is  it  possible  for  the  shuttle  to  make  an  abort  landing  at  an  East  coast  airport? 
This  thesis  will  explore  the  possibilities  of  safely  landing  the  Orbiter  at  an  East  Coast  airport  by  first 
creating  a  computer  model  of  a  shuttle  laimch  from  5 1.6  degrees  and  57  degrees  then  designing  another 
computer  model  using  Hypersonic  Orbiter  reentry  data  collected  from  NASA  to  mimic  the  dynamics  of  the 
shuttle  making  an  abort  landing.  The  Launch  model  will  give  the  Abort  model  initial  state  vectors  from  4 
minutes  20  seconds  until  just  prior  to  main  engine  cutoff.  Several  control  methods  will  be  used  to  land  the 
Orbiter.  The  first  method  will  use  constant  pitch  and  roll  until  the  Orbiter  reaches  a  preset  altitude.  The 
second  will  hold  pitch  and  roll  constant  but  change  pitch  to  30  degrees  and  roll  to  0  degrees  when  the 
desired  heading  angle  is  achieved  to  intercept  a  target  airport.  Both  methods  will  use  the  same  abort  model 
with  minor  changes  in  output  and  the  afore  mentioned  control  commands.  The  first  method  will  map  out 
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possible  areas  of  intercept  over  the  map  of  the  East  Coast  as  well  as  final  state  vectors  giving  longitude, 
latitude,  speed,  and  altitude.  The  location  with  respect  to  airports  of  sufficient  runway  lengths  and  low 
enough  speeds  with  respect  to  nominal  re-entries  will  prove  the  feasibility  of  this  idea  as  well  as  give  an 
idea  as  to  which  airports  the  second  method  will  use  as  target  landing  sites.  The  second  method  will  map 
out  all  abort  descents  throughout  the  launch  trajectoiy  to  the  target  airports  selected  from  the  first  method. 
This  will  give  clear  times  throughout  the  launch  when  certain  airports  become  abort  landing  sites. 
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2.  Methodology 


2.1  Introduction 

This  chapter  will  explain,  in  detail,  how  the  problem  of  examining  possible  abort  landings  at  east- 
coast  airports  from  launches  of  51.6  and  57  degree  inclinations  out  of  Cape  Canaveral  is  approached  and 
resolved  using  computer  modeling  in  FORTRAN,  The  first  section  explains  how  the  launch  of  the  Shuttle 
was  computer  modeled  as  a  gravity-turn  trajectory.  This  program  provides  state  vectors  at  every  point 
throughout  the  launch,  which  are  used  as  initial  points  in  the  abort  landing  program.  The  second  section 
discusses  the  abort  landing  computer  model  that  uses  hypersonic  equations,  as  well  as  NASA  data  from 
Shuttle  re-entries  to  take  the  orbiter  without  external  tank  and  bring  it  to  a  nominal  landing  altitude  using 
two  methods  of  controlling  the  descent  using  pitch  and  roll  without  thrust.  For  both  methods  pitch  is 
measured  as  positive  nose  up  and  roll  is  measured  positive  roll  left.  The  first  method  uses  constant  pitch 
and  roll  angles  for  the  entire  descent  while  the  second  method  only  holds  pitch  and  roll  constant  until  a 
certain  heading  angle  is  reached  then  roll  becomes  0  degrees  and  pitch  becomes  30  degrees.  Output  data 
for  both  methods  are  mapped  against  airport  locations  which  have  sufficient  runway  lengths  to 
accommodate  the  Space  Shuttle.  These  two  control  methods  will  give  a  clear  answer  to  the  possibility  of 
an  east  coast  abort  mode. 

2.2  Space  Shuttle  Launch  Model 

2.2.1  Assumptions.  When  attempting  to  model  reality,  the  first  step  is  to  minimize  the  number  of 
variables  so  that  the  problem  becomes  workable.  However,  eliminating  too  many  variables  will  trivialize 
the  problem  and  the  model  will  be  nothing  close  to  reality.  For  the  launch  program  a  flat  non-rotating  Earth 
was  assumed.  This  sounds  like  a  huge  departure  from  reality,  but  when  looking  at  the  problem  from  the 
point  of  view  of  what  is  needed,  i.e.  the  state  vector,  the  only  states  which  involve  a  round  Earth  are 
longitude  and  latitude.  These  states  can  be  derived  from  the  downrange  distance  using  spherical 
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trigonometry  and  the  acceleration  terms  in  the  equations  of  motion  are  modified  to  include  the  apparent 
“centrifugal  force.”  So,  in  effect,  a  spherical  problem  is  made  flat. 

The  second  assumption  is  that  the  velocity,  thrust,  and  drag  vector  lie  along  the  same  line.  This 
simply  means  the  Shuttle  always  travels  in  the  direction  its  nose  is  pointing.  Furthermore;  the  tlirust  is 
assumed  to  be  constant  and  perfectly  efficient.  This  eliminates  both  the  possibilities  of  throttling  the 
engines  and  any  thrust  loss  due  to  change  in  altitude.  The  coefficient  of  drag  is  assumed  to  be  1 .0  for  all 
components  of  the  space  shuttle  system.  Of  course,  each  component  will  have  different  values  of  drag,  due 
to  the  differences  of  fi-ontal  surface  area. 

Two  more  assumptions  are  made  at  the  beginning  of  the  launch  since  the  initial  state  of  the 
program  starts  when  the  Shuttle  is  already  100  meters  off  the  ground.  Time  and  mass  start  at  0  and  MO 
respectively,  meaning  that  it  took  no  time  to  fly  the  first  1 00  meters  and  no  mass  was  lost  getting  there. 
Since  this  is  a  multi-stage  launch,  it  is  also  assumed  that  solid  rocket  booster  and  external  tank  separation 
happen  at  exactly  120  seconds  and  480  seconds,  respectively. 

The  final  assumption  is  that  the  atmosphere  behaves  exactly  like  the  Regan  and  Anandarskarian 
Atmosphere  Model.  Figure  2-1  describes  the  density  profile  of  this  model  which  shows  the  density  to  be 
substantially  less  above  50  km  altitude  (9:  Appendix  A) 
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2,2.2  Statistics  Used for  Launch .  Table  2-1  lists  Shuttle  mass  and  performance  data,  as  well  as 
constants  necessaiy  for  the  computer  program. 


Table  2-1.  Shuttle  Statistics  ( 7:490;  6:1 


Shuttle  Mass 

Orbiter  -f  Payload 

108,864  Kg 

2  Solid  Rocket  Boosters(SRBs) 

164,000  Kg 

External  Tank 

32,000  Kg 

Fuel 

1,736,336  Kg 

Total  Launch  Mass  ( MO) 

2,041,200  Kg 

Shuttle  Petformance 

Thrust  at  Launch 

30,300kN 

Thrust  at  Second  Stage 

6,300kN  (at  95%  6,000kN) 

Mass  Flow  Rate  at  Launch 

9,906  Kg/s 

Mass  Flow  Rate  at  Second  Stage 

1,496  Kg/s 

Constants 

Radius  of  the  Earth  ( R  ) 

6,378.135  Km 

Gravitational  Parameter  of  the  Earth  (p ) 

3.98601E14mVs^ 

Sea  Level  Pressure  (  PO  ) 

101,325  N/m^ 

Front  Surface  Area  of  the  Orbiter 

29.235  (Estimated) 

Front  Surface  Area  of  External  Tank 

50.265  m^ 

Front  Surface  Area  of  Both  SRBs 

21.621  m^ 

Coefficient  of  Drag  (  Cd  ) 

1.0  (Assumed) 

Longitude  of  Cape  Canaveral 

-80.60  degrees 

Latitude  of  Cape  Canaveral 

28.46  degrees 

Initial  Launch  Conditions 

Down  Range  Distance  ( X ) 

Om 

Altitude  (H) 

100  m 

Velocity  (V) 

32.1  m/s 

Flight  Path  Angle  (y ) 

1.56209405  radians 

Mass  (  M) 

2,041,200  Kg 
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2.2.3  Gravity-Turn  Trajectory.  The  dynamics  of  the  Space  Shuttle  launch  model  were  chosen  to 
match  the  equations  governing  a  gravity  turn  trajectory.  This  type  of  launch  aligns  the  thrust  and  velocity 
vectors,  as  mentioned  in  the  assumptions,  as  well  as  the  drag  vector.  This  is  done  primarily  for  launch 
vehicles  that  cannot  sustain  large  aerodynamic  loads  in  the  transverse  direction,  but  can  carry  large  loads  in 
the  axial  direction  (13:207). 


Figure  2-2.  Free  Body  Diagram  of  Launch  Vehicle  (13:208) 


Figure  2-2  shows  the  free  body  diagram  of  the  launch  vehicle  with  all  forces  acting  on  it.  The  equations  of 
motion  for  this  type  of  launch  trajectory  are  as  follows: 


dt 


=  V  cos;" 


(2.1) 


dH 


—  ^Vsinr 


(2.2) 
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(2.3) 


dV  ^  r.  r  rnX\  . 
^-J-T-D-irng-—)smy 


u‘‘y  <  ^ 


(2.4) 


where 

X-  Downrange  Distance  (m) 

H- Altitude  (m) 

V-  Speed  (m/s) 

y=^  Flight  Path  Angle  (radians) 

M=  Mass  (kg) 

7^  Thrust  (Newton’s) 

D  =  Drag  (Newton’s) 

R  =  Radius  of  the  Earth  (m) 
t  =  Time  (s) 

The  only  noticeable  irregularity  in  these  equations  are  the  bracketed  terms  in  equations  2.3  and  2.4.  These 
terms  include  apparent  “centrifl^al  force”  so  that  flight  over  a  spherical  Earth  can  be  resolved  to  one  over  a 
flat  Earth  (13:  208).  These  equations  must  be  modified  further  to  accurately  model  an  actual  launch  with 
statements  pertaining  to  the  change  in  mass,  thrust,  and  drag  with  time.  Both  mass  flow  and  thrust  were 
obtained  from  Space  Shuttle  data  and  the  drag,  which  is  discussed  later  in  this  chapter,  was  estimated. 

The  first  step  in  numerically  integrating  these  equations  is  to  find  the  initial  conditions.  The  initial 
launch  can’t  start  with  a  flight  path  angle  of  90  degrees.  Otherwise,  by  equation  2.4,  y  would  never  change. 
This  is,  of  course,  not  desired  so  a  slight  “nudge”  is  necessary  to  cause  y  to  change  and  by  the  negative  sign 
in  equation  2.4  it  will  decrease,  causing  the  vehicle  to  fall  over.  If  this  is  done  on  the  pad,  the  booster 
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would  be  lying  on  the  ground  in  a  matter  of  seconds.  To  prevent  this  from  happening,  the  pitch  over  or 


“nudge”  is  performed  at  100  meters  altitude.  This  is  chosen  because  it  is  best  to  do  the  initial  pitch  over  at 
low  speed,  since  the  vehicle  will  briefly  have  an  angle  of  attack  during  the  maneuver,  and  high  speeds 
would  inflict  greater  aerodynamic  loads. 

After  the  initial  altitude  is  chosen,  the  velocity  can  be  determined  by  integrating  the  launch  in  a 
straight  line  from  0  to  100  meters.  Knowing  the  thrust,  mass,  and  drag  on  the  space  shuttle,  the  initial 
velocity  was  calculated  to  be  32.1  m/s.  As  mentioned  in  the  assumptions,  mass  was  assumed  not  to  change 
over  this  period.  Now  that  all  the  initial  conditions,  except  for  y  have  been  found,  the  only  thing  left  to  do 
is  set  final  conditions  and  iterate  different  initial  values  of  y  until  those  final  conditions  are  met. 

Final  conditions  were  established  at  each  stage  of  the  Shuttle’s  launch  profile  by  data  from 
NASA.  At  the  first  stage,  altitude  was  given  to  be  50  km  and  velocity  had  to  be  1.4  km/s  (10:  4).  At  main 
engine  cut-off,  y  had  to  be  0,  while  altitude  and  velocity  had  to  be  145  km  and  7.82  km/s,  respectively.  The 
velocity  was  calculated  from  the  altitude  of  145  km  using  the  equation  for  bum-out  velocity  of  a  circular 
orbit  given  by 


Vbumout  —  ( 


R  +  Hmeco 


1/2 


(2.5) 


where 

jii  =  Gravitational  Parameter  for  the  Earth  ( m  W) 

R  -  Radius  of  the  Earth  (m) 

^meco  =  Altitude  at  Main  Engine  Cut-Off  (m) 

These  conditions  were  used  as  a  reference  to  model  file  launch  so  that  optimally  the  launch  should 
look  like  Figure  2-3(13: 210). 
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Figure  2-3.  Gravity-Tum  Trajectory 


However,  when  the  equations  of  motion  were  numerically  integrated  an  angle  for  y  could  not  be  found  to 
match  both  staging  conditions.  This  is  most  likely  do  to  the  fact  that  in  reality  the  Shuttle  throttles  it’s 
engines  and  since  this  model  uses  constant  thrust  a  slight  modification  must  be  made.  A  control  input  for  y 
is  used  at  various  points  to  achieve  both  conditions  so  the  actual  computer  model’s  launch  looks  like  Figure 


2-4. 
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Figure  2-4.  Modified  Gravity-Turn  Trajectory 


This  model  doesn’t  exactly  match  the  optimum  trajectory  but  provides  a  close  enough  approximation  to 
give  reasonable  state  vectors  necessary  for  the  abort  program.  The  value  for  the  initial  y  was  found  to  be 
1.56209405  radians  with  an  increase  of  .1 1  radians  at  210  seconds,  .12  radians  at  280  seconds,  .10  radians 
at  340  seconds,  and  .02  radians  at  420  seconds. 

2.2.4  Drag.  Before  numerical  integration  is  possible,  it  is  first  necessary  to  define  the  drag  term 
in  equation  2.3.  The  drag  law  states  that 
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D  =  ^CdApV^ 


(2.6) 


where 


Q)  =  Coefficient  of  drag  (assumed  to  be  1 .0) 

A  =  Frontal  surface  area  (m^  ) 

V-  Speed  (m/s) 

Po  =  Atmosphere  density  (kg/m^) 

The  frontal  surface  areas  of  the  SRBs  and  external  tank  were  calculated  from  the  equation 


A  = 


(2.7) 


where 


d  -  Diameter  (m) 

The  frontal  area  of  the  Orbiter  could  only  be  estimated  by  adding  the  area  of  both  the  fijselage  and  the 
wings.  The  frontal  area  of  the  wings  were  estimated  by  multiplying  the  wingspan,  24m,  by  an  assumed 
value  for  the  wing  thickness,  0.4  m. 

The  atmosphere  density  was  found  using  a  FORTRAN  subroutine  called  Atm,  which,  given 
altitude  and  ground  level  pressure,  outputs  density  at  altitude.  This  subroutine  uses  the  atmosphere  model 
found  in  the  appendix  of  Regan  and  Ananderskarian,  as  mentioned  in  section  2.2.1. 

The  other  parameter,  coefficient  of  drag,  is  also  mentioned  in  2.2. 1 .  Although  Cd  is  known  to  be  a 
fimction  of  mach  number,  Reynolds  number,  and  aerodynamic  configuration,  it  is  assumed  as  constant. 

The  values  of  Cp  may  range  from  0. 12  to  2.0,  depending  on  the  aforementioned  conditions,  however  1 .0 
was  selected  as  a  mid-range  value  as  well  as  for  simplicity  (1 :55). 

2.2.5  Numerical  Integration.  The  equations  of  motion  and  initial  conditions  for  the  gravity-turn 
are  now  known.  What  is  needed  now  is  to  first  put  equations  2. 1-2.4  into  the  form 

X  =  f{x,t)  (2.8) 

which  is  done  simply  by  dividing  equation  2.3  by  the  mass  and  equation  2.4  by  both  mass  and  velocity. 

The  next  step  is  to  ensure  that  the  state  vector  elements  have  close  to  the  same  order  of  magnitude.  SI  units 
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were  chosen  for  both  their  ease  in  conversions  and  proper  orders  of  magnitude.  The  last  step  is  then  to  find 
a  numerical  integrator  which  is  best  suited  for  the  particular  problem. 

Numerical  integrators  fall  into  several  classes:  extrapolators,  predictors,  and  predictor-correctors. 
Extrapulators  and  predictors  step  their  way  into  the  future  using  data  from  the  current  instant  and  the 
immediate  past  (12:1 19).  Predictor-correctors  add  a  corrector  algorithm  which  runs  a  higher  order 
polynomial  through  the  previous  data  points  to  obtain  a  better  value  of  the  new  state  vector  (12:123).  From 
these  choices,  a  fourth  order  predictor-corrector  subroutine  in  FORTRAN  called  Haming  was  selected 
(Appendix  A).  Although  Haming  is  over  fifteen  years  old,  it  still  performs  better  than  ODE  on  the  VAX  or 
CYBER  by  a  substantial  margin  (12:120).  Since  this  integrator  is  fourth  order,  the  only  drawback  is  that 
the  initial  state  provides  only  one  point  and  an  initialization  method  is  needed  to  get  the  other  three  before 
they  can  be  called  into  operation. 

In  addition,  like  most  integrators,  Haming  prefers  smooth  functions  and  since  this  is  a  multi-stage 
launch  where  mass  and  thrust  change  instantaneously  at  various  stages,  Haming  must  be  re-initialized  at 
each  stage.  To  accomplish  this,  it  was  necessaiy  to  write  three  different  subroutines  called  Rhsl,  Rhs2, 
Rhs3,  which  all  include  equations  2. 1-2.4  but  have  different  thrust  and  mass  values  corresponding  to  stage. 
These  three  subroutines  are  then  integrated  using  integrators  Hamingl,  Haming2,  Hammg3  and  added 
together.  This  provides  a  smooth  integration  over  a  function  with  several  jumps. 

2.2.6  Spherical  Trigonometry.  The  original  state  vector,  [X,  H,  V,  y ,  M],  having  gone  into  the 
program  and  been  numerically  integrated  is  output  as  that  state  at  any  time  along  its  flight.  However,  the 
desired  output  is  to  have  location  measured  in  longitude  and  latitude  on  the  globe  rather  than  downrange 
distance,  X.  We  also  want  to  know  the  heading  angle  at  every  instant  after  launch.  This  step,  in  effect, 
turns  the  2D  launch  into  three  dimensions.  Latitude,  longitude,  and  heading  angle  can  be  computed  using 
the  original  outputs  and  few  added  initial  states. 

The  first  step  to  modifying  these  outputs  is  to  calculate  the  initial  heading  angle,  ij/o ,  given  the 
inclination  of  the  launch  and  the  longitude  and  latitude  of  the  launch  site.  Figure  2-5  displays  a  right 
triangle  made  from  three  great  circles  (8:20). 
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The  two  vertices  that  are  not  90  degrees  represent  the  equator  and  the  launch  site,  and  the  inscribed  angles, 
A  and  B,  represent  the  launch  inclination  and  azimuth  angle  respectively.  The  two  sides,  a  and  b,  are  the 
latitude  of  the  launch  site  and  the  longitudinal  distance  from  the  equatorial  crossing.  From  Napier’s  Rule,  it 
is  foxmd  that 


COS  A  =  cosa  sin  B  (2.7a) 

or  rearrange 


B  =  sin  ' 


cosv4 

cosa 


(2.7b) 


The  initial  heading  angle  can  then  be  found  by  subtracting  the  azimuth  from  90  degrees. 

The  next  step  is  to  find  the  longitude  and  latitude  throughout  the  flight,  given  the  conditions 
already  stated  and  the  downrange  distance. 
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Figure  2-6  shows  yet  another  spherical  triangle  with  the  three  vertices  as  the  launch  site,  Orbiter,  and  North 
Pole.  Side  c  is  the  downrange  distance  converted  to  radians  by  dividing  by  the  radius  of  the  earth  making  it 
a  great  circle.  The  other  known  values  are  a ,  which  is  the  difference  of  90  degrees  and  the  latitude  of  the 
launch  site,  and  B  which  is  the  launch  azimuth  found  in  equation  2.7b.  The  following  equations  are  used  to 
solve  for  the  unknowns 


b  =  cos  '  [cosa  cos  c  +  sin  a  sin  c  cos  5]  (2.8) 


C  =  cos  ’ 


cosc-cosbcosa 

sin^sina 


(2.9) 


where 


b  =  Difference  of  90  Degrees  and  the  Latitude  of  the  Orbiter 
C=  Difference  in  the  Longitude  of  the  Launch  Site  and  the  Longitude  of  the  Orbiter 
Then  all  that  remains  is  to  solve  both  b  and  C  for  the  latitude  and  longitude  of  the  Orbiter. 

iaUter=90-b  (2.10) 
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Loflgorbiter  ^^^Ssite  ”  ^ 


(2.11) 


where 

Latsite  =  Latitude  of  launch  site 
Longsite  =  Longitude  of  launch  site 

The  heading  angle  can  then  be  evaluated  at  every  point  in  the  flight  by  subtracting  equation  2.7b  from  90 
degrees  with  ‘a’  as  the  latitude  of  the  launch  vehicle,  rather  than  that  of  the  launch  site. 

2.2.7  Flow  Chart  for  Gravturn  Program,  Figure  2-7  shows  how  each  subroutine  mentioned  in  the 
previous  sections  of  this  chapter  interface  with  each  other  to  provide  the  necessary  output. 


FLOW  CHART  FOR  GRAVTURN  PROGRAM 


Figure  2-7.  Flow  Chart  of  Gravturn  Program 


The  inputs  are  the  angle  of  inclination,  along  with  the  initial  state  vector:  downrange  distance,  altitude, 
velocity,  flight  path  angle,  and  mass.  The  output  is  the  new  state  vector  at  any  time  during  the  flight: 
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altitude,  longitude,  latitude,  velocity,  flight  path  angle,  and  heading  angle.  The  main  program  and  all 
subroutines  can  be  found  in  Appendix  A  of  this  thesis. 

2.3  Abort  Landing  Model 

Many  different  circumstances  may  cause  the  Shuttle  to  choose  an  abort  mode.  These  modes  vary 
depending  on  which  systems  have  failed  and  if  time  to  land  is  critical.  This  thesis  cannot  explore  every 
possible  combination  of  working  and  malfunctioning  systems,  but  it  will  attempt  to  resolve  the  feasibility  of 
aborting  without  any  working  propulsion  system  given  the  flight  control  system  still  operates.  The  abort 
program  will  take  the  Shuttle  at  any  point  in  it’s  launch  after  SRB  separation,  drop  the  external  tank,  and 
glide  in  using  only  pitch  and  roll  to  control  the  un-powered  descent.  This  gives  a  clear  idea  as  to  which 
times  throughout  the  launch,  if  any,  it  will  be  possible  to  land  at  an  east  coast  airport. 

2,3.1  Assumptions,  For  this  computer  model,  the  thrust  is  assumed  to  be  zero.  Since  the  external 
tank  is  dropped,  the  mass  will  be  only  that  of  the  Orbiter  plus  payload.  The  propulsion  system  is  non¬ 
operative  and  cannot  expend  fuel,  meaning  the  mass  is  constant.  The  control  system  on  board  is  operative, 
which  leads  to  the  largest  assumption,  an  aerodynamic  model  of  the  Orbiter. 

The  Rarefied  -  Flow  Shuttle  Aerodynamics  Flight  Model  will  be  assumed  to  provide  Cl  and  Cd 
values  for  varying  altitudes,  speeds,  and  vehicle  attitudes.  This  model  is  taken  from  NASA  data  collected 
on  twelve  Orbiter  re-entries  and  will  be  discussed  later  in  this  chapter.  Although  the  field  of  rarefied-flow 
is  highly  theoretical,  it  is  based  on  actual  data  and  should  serve  as  the  best  assumption  available  for 
modeling  winged  vehicles  through  atmosphere  at  hypersonic  speeds. 

Further  assumptions  are  made  to  use  equations  based  on  flight  over  a  spherical  planet.  One  is  that 
the  planet  is  rotating  at  such  a  rate  that  the  coriolis  acceleration  term,  2coV,  is  retained  but  the  higher  order 
terms  such  as  o^r  are  neglected.  The  second  is  that  the  Orbiter  is  treated  as  a  point  mass.  This  only  means 
the  forces  acting  on  the  Orbiter  are  of  primary  interest  whereas  the  moments  can  be  looked  at  in  a  more  in 
depth  study  if  east  coast  aborts  are  indeed  considered  feasible. 
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2.3.2  Statistics  Used  For  Shuttle  Abort  Model.  Table  2-2  shows  the  values  used  to  represent  the 


performance  and  specifications  of  the  Shuttle  for  the  Abort  Program. 


Table  2-2.  Data  Used  for  Abort  Program  ( 3:550-555;  2:429;  13:316 ) 


Shuttle  Data 

Orbiter  &  Payload  Mass  (M) 

108,864  Kg 

Vehicle  Reference  Area  (S) 

249.9 

Mean  Aerodynamic  Chord  ( MAC  ) 

12.058  m 

Constants 

Radius  of  the  Earth  ( R  ) 

6378.135  Km 

Gravitational  Parameter  of  the  Earth*  (  p) 

3.98601x10’'*  mV 

Rotation  Rate  of  the  Earth  (co) 

7.2921 15856x10-®  radians/s 

Sea  Level  Pressure  (Po  ) 

101325  N/m^ 

Final  Conditions  for  Nominal  Landing 

Altitude 

25,000  m 

Velocity 

800  m/s 

2.3.3  Flight  Over  A  Spherical  Planet.  This  section  will  derive  and  explain  the  equations  of 
motion  (EOMs)  that  were  used  in  the  abort  landing  program.  As  mentioned  in  the  assumptions,  these 
equations  govern  the  motion  of  a  point  mass  with  zero  thrust  flying  inside  an  atmosphere  of  a  spherical 
planet.  They  describe  the  time  rate  of  change  of  the  radius,  longitude,  latitude,  velocity,  flight  path  angle, 
and  heading  angle. 
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Figure  2-8  shows  a  planet  fixed  system  XYZ  and  another  system  xyz  which  is  rotating  with  respect 


to  it. 


Figure  2-8.  Coordinate  Systems  (1 1:22) 


The  planet  fixed  system  is  rotating  with  the  angular  velocity  of  the  earth,  o)^  assumed  constant  and  directed 
along  the  Z  -  Axis  (1 1:21).  The  vectors  r  and  v  represent  the  position  and  velocity  of  the  point  mass  M, 
The  angles  Band  ^  represent  the  longitude  (  measured  from  the  X  Axis,  in  the  equatorial  plain,  and 
positively  eastward)  and  latitude  (measured  from  the  equatorial  plane,  along  a  meridian,  and  positively 
northward)  respectively.  As  in  the  gravity-turn  trajectory  y  and  y/  are  the  flight  path  angle  and  heading 
angle  respectively.  The  flight  path  angle  is  measured  between  the  local  horizontal  and  the  velocity  vector 
while  the  heading  angle  is  between  the  local  parallel  of  latitude  and  the  projection  of  v  on  the  horizontal 
plane.  The  angle  o’ or  roll  angle  is  measured  from  the  direction  of  lift  from  the  Orbiter  and  the  plane  made 
by  the  position  and  velocity  vector  (11:  23-24). 

If  i ,  j  ,  and  k  are  unit  vectors  along  the  xyz  reference  frame,  the  position,  velocity,  and  angular 
velocity  of  the  earth  can  be  written  in  component  notation  as 
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r  =  ri 


(2.12) 


v  =  vsin;'  /  +  vcos;' cos^^  7  +  vcos;' k  (2.13) 

m-co  s\n<f>  I  +(ocos^  k  (2.14) 

where 

r  =  Magnitude  of  the  Radius  (m) 

V  =  Magnitude  of  the  Velocity  (m/s) 

6>  =  Angular  Velocity  of  the  Earth  (radians/s) 

Y  =  Flight  Path  Angle  (radians) 

1!/=  Heading  Angle  (radians) 

Latitude  (radians) 


It  follows  that  acceleration  can  be  found  using  the  transport  theorem 


dv  dv 

inertial  —  local  +  X  V  CO  x(jCO  X  v)  (2.15) 


where 


dv 

local  =  Time  Rate  of  Change  of  Velocity  wrt  the  Local  Reference  Frame 
at 


Substituting  Eq  2.15  into  Newton’s  second  law  results  in 


M  local  —  F-2Ma  X  v  —  M(o  x  (co  x  r) 
dt 


where 


M  =  mass  (kg) 

F  =  Total  Force  Vector  (Newton’s) 


(2.16) 
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The  total  force  is  represented  in  the  free  body  diagram  of  figure  2-9.  The  force  of  gravity,  mg,  is  in  the 
negative  x  direction  while  the  force  of  drag  is  in  the  negative  direction  of  the  velocity  vector.  All  that 
remains  to  have  a  complete  knowledge  of  all  the  forces  involved  is  to  determine  the  direction  of  lift  in  the 
xyz  coordinate  frame. 


Diagram  of  Forces 


Figure  2-9.  Force  Diagram 

Let  XiyiZj  be  an  intermediate  coordinate  system  rotated  an  angle  a  firom  a  body  fixed  frame  as 
shown  in  figure  2-10. 
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The  system  XiyiZi  is  deduced  by  a  rotation  \\f  in  the  horizontal  plane  then  a  rotation  y  in  the  vertical  plane 
of  the  xyz  system  (11 :24).  The  transformation  matrix  equation  for  this  sequence  of  rotations  is 


X 

cosy 

sin;^ 

0 

Xi 

y 

- 

-siny  cosy/ 

cos/ 

z 

-sinz  siny/ 

cosy 

cosy/ 

Zi 

which  resolves  the  direction  of  lift  in  the  xyz  system  to 

L  =  Lcosa  cosy  i -(Looser  sin;^  cos^^  +  Zsincr  smyr)] 
-(Looser  siny  siny/  +  Lsiner  oosy/)k 


(2.18) 


The  total  force,  position,  velocity,  and  angular  velocity  of  the  earth  are  now  known  in  terms  of  the  unit 
vectors,  i,  j,  k  in  the  xyz  frame.  In  order  to  take  the  derivatives  of  the  vectors  with  respect  to  the  XYZ 
fi*ame,  the  derivatives  of  the  unit  vectors  i,  j,  k  must  be  determined. 

In  much  the  same  fashion  as  the  earth  fixed  referenced  frame  XYZ  was  evaluated  in  an  inertial 
frame  by  means  of  the  transport  theorem  via  the  rotation  rate  of  the  earth,  (d  ,  the  xyz  frame  must  be 
evaluated  via  its  rotation  rate,  relative  to  the  inertialy  fixed  planet.  The  system  xyz  is  acquired  from  the 
XYZ  system  by  a  rotation  9  about  the  positive  z  axis  then  by  a  rotation  ^  about  the  negative  y  axis  (11 :24). 
This  yields  the  angular  velocity 

Q  =  (sm^— )/  -  (— )y  +  {cos<f)—)k  (2.19) 

The  time  derivatives  can  then  be  obtained  using  Poisson  formulas  (1 1:20) 
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-=n><,=(cos«>-)y+(-)t 


—  =  Q  X  y  =  -(cos(^— )/  +  (sin^>— )^ 


(2.20) 


-=nxi  =  -(-),-(.m^-); 


The  next  step  is  to  take  the  time  derivatives  of  equations  2.12  and  2.13,  which  for  position  yield 

dr  dr  ^ 

inertial  -  +  (''  COS^— )7  +  (t  (2-21) 


which  can  be  broken  up  by  component  unit  vectors  and  set  equal  to  equation  2.13  yielding 


dr 

-  =  vsmr 


d9  vcosycos^ 
dt  rcos^ 


(2.22) 

(2.23) 


d<^  vcosysmy/ 
dt  r 


(2.24) 


Eqriations  2.22, 2.23,  and  2.24  represent  the  first  three  equations  of  motion  used  in  the  abort  program.  The 


second  three  equations  of  motion  are  obtained  from  taking  the  derivative  of  equation  2. 13  and  substituting 


dr  dO  d^ 

the  first  three  EOMs  for  “ ,  ,  and  .  This  resxilt  along  with  the  values  for  total  force  and  co  are 

dt  dt  dt 
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dv  dy 

placed  into  Eq  2.16,  the  force  equation.  This  equation  is  then  rearranged  to  solve  for  —r,  — 

dt  dt 


,  and 


dt 


After  neglecting  terms  which  have  (o\,  the  final  three  EOMs  are 

dv  -D 


dy  _  Lcosa 
dt  Mv 


gcosy  vcosy 

- + - +  26?C0S®C0SI/ 

V  r 


(2.26) 


dy/  Lsina  vcosy  cosy/tm^ 

- ~ - +  2G>(tanrcos^sinMr-sin^)  (2.27) 

dt  Mvcosy  r 


where 


D  =  drag  (Newton’s) 
L  =  lift  (Newton’s) 
M  =  mass  (kg) 
g  =  gravity  (m/s^) 


These  six  equations  were  placed  in  a  subroutine  called  Dynam  where  they  can  be  numerically  integrated  by 
the  Haming  subroutine  used  in  Gravtum  once  the  values  for  Drag  and  Lift  are  determined. 

2.3.4  Shuttle  Aerodynamics  Flight  Model  Equations  2.22-2.27  describe  the  dynamics  for  the 


abort  model  and  are  known  but  can  not  be  numerically  integrated.  Equations  2.25, 2.26,  and  2.27  have  lift 
and  drag  terms  whose  equations  are 


(2.29) 


1  , 
D  =  —Cdp  v^S 


where 

Cl  =  coefficient  of  lift 
Cr>  =  coefficient  of  drag 
p  =  density  (kg/m^) 

V  ==  velocity  (m/s) 

S  =  reference  area  (m^) 

Density  is  found  using  the  same  earth  atmosphere  model  of  figure  2-2.  However,  the  coefficients  of  lift  and 
drag  can’t  be  assumed  constant;  therefore  a  model  must  be  found  to  accurately  estimate  the  lift  and  drag 
coefficients  the  Orbiter  will  experience  at  various  altitudes  and  speeds.  The  obvious  choice  is  the  Shuttle 
rarefied  -  flow  aerodynamics  flight  model  that  was  derived  from  a  High  Resolution  Accelerometer  Package, 
HiRAP,  on  twelve  Shuttle  re-entries  spanning  a  period  of  ten  years  (6:550). 

Rarefied  -  flow  is  the  transition  region  between  fi'ee  molecular  flow  and  hypersonic  continuum, 
which  is  roughly  160-60  km  in  altitude.  These  regions  are  determined  based  on  the  Knudsen  number,  Kn, 
Of  ratio  of  the  mean  free  path  to  the  mean  aerodynamic  chord,  MAC  (6:  552).  The  regions  are  clearly 
separated  in  Table  2-3. 


Table  2-3.  Atmospheric  Regions  Based  on  Knudsen  Number  (3:550) 


Kn 

Region 

Hypersonic  Continuum 

10'^<Kn<10 

Rarefied-Flow 

>10 

Free  Molecular  Flow 
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Since  the  majority  of  the  time  the  Orbiter  will  be  in  the  rarefied  -  flow  region  it  will  be  necessary  to  use  the 
empirical  equations  which  govern  it.  The  equations  are  split  into  three  segments:  Hypersonic  continuum. 
Free  molecular  flow,  and  a  bridging  formula.  All  three  segments  calculate  Cn  and  Ca,  which  are 
aerodynamic  coefficients  in  the  normal  and  axial  direction,  as  functions  of  pitch  angle,  a.  The  reason  for 
using  normal  and  axial  rather  than  lift  and  drag  is  primarily  due  to  the  placement  of  the  HiRAP  on  board  the 
Orbiter,  This  only  adds  one  further  step  which  will  convert  Cn  and  Ca  to  Cl  and  Cq. 

The  hypersonic  continuum  equations  for  the  normd  and  axial  coefficients  as  a  function  of  pitch 
are 


Cnc  =  -9.25704  X 10"^ +  5.23808  x  1 0'^  a  -  0.839782 


(2.30) 


C.<c  =  5.86689  X  10"’a^  - 6.72027  x  1 0“^ cr"  +3.32044  xlO'^ a -0.00863 14  (2.31) 

these  equations  are  curve  fit  to  match  Shuttle  data  for  an  envelope  of  35<a<45  degrees  (3:552). 

The  free  molecular  flow  equations  for  Cn  and  Ca  as  a  function  of  a  are 


Cnf  =  -7.16528  X  lO-'^a^  +  9.66197  x  lO'^^  +  9.18422  x  lO'^a  +  li8739  x  10'^  (2.32) 

= -1.171 17  xlO'^cr^  +  5.92205  xlO^'a^  +  0.0164864a +  0.751 105  (2.33) 

These  equations  were  curve  fit  to  match  pitch  angles  between  0  and  60  degrees  (3:553).  All  that  remains  to 
compute  Cn  and  Ca  for  the  Shuttle  is  the  bridging  formula  to  combine  hypersonic  continuum  with  free 
molecular  flow. 

The  bridging  formula  depends  on  the  value  of  Kn  as  shown  in  the  following  equations 
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Cn  =  exp[-  0.29981(1.3849  -  logio^n)*^*"*] 

if  log  loATn  <  1,3849  (2.34) 

otherwise  Cn  =  1.0 


Ca  =  exp[-  02262(1.2042  -  log  io^«)  * 

if  log  ioi&7  <  1.2042  (2.35) 

otherwise  Ca  =  1.0 


Now  the  results  from  equations  2.30  -  2.35  are  used  to  calculate  the  final  normal  and  axial  aerodynamic 
coefficients  of  the  Orbiter 

Cn  =  Cnc  +  (Cnf  —  Cnc)Cn  (2.36) 

Ca  =  Cac  +  (Cap  -  Cac)Ca  (2.37) 

The  normal  and  axial  coefficients  can  now  be  converted  to  Cl  and  Cn-  Figure  2-1 1  shows  a  vector 
representation  of  the  alignment  of  these  coefficients  on  a  vehicle  with  an  angle  of  attack. 
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Diagram  of  Aerodynamic  Coefficients 


Figure  2-11.  Vector  Components  of  Aerodynamic  Coefficients 


The  dashed  components  make  up  Cl  and  Cd 


Cl  =  -Ca  sin  a +  Cn  cos  a  (2.38) 

Cd  =  Ca  cos  a +  Cn  sin  a  (2.39) 


The  equations  from  2.30  to  2.39  were  placed  in  a  subroutine  called  Aero  which  returns  Cl  and  Cd  given 
pitch  and  Kn.  This  in  turn  provides  the  EOMs  placed  in  the  subroutine  Dynam  with  lift  and  drag  so  that 
they  may  now  be  numerically  integrated  much  the  same  way  the  gravity  turn  trajectory  equations  were. 
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2.3.5  Flow  Chart.  Figure  2-12  shows  how  the  abort  program  interacts  with  it’s  subroutines. 


Flow  Chart  of  Abort  Program 


Constants 


Figure  2-12.  Flow  Chart  Of  Abort  Program 


The  inputs  to  the  program  are  the  state  vector  produced  by  the  Gravtum  Program,  the  constants  necessary  to 
the  reentry  equations,  and  varying  values  of  pitch  and  roll.  The  program  then  calls  Haming  which  in  turn 
calls  D5mam.  Dynam  contains  the  six  EOMs  that  describe  flight  over  a  spherical  planet  which  need  density 
and  the  coefficients  of  lift  and  drag  that  it  receives  from  Atm  and  Aero  respectively.  The  output  of  the 
program  provides  a  final  longitude,  latitude,  altitude,  and  velocity  which  are  mapped  against  existing 
airports  to  determine  if  east  coast  aborts  are  possible. 
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3.  Data  Analysis  and  Results  for  Constant  Control  Descent 


3,1  Introduction 

This  chapter  presents  the  data  from  the  Abortl  program  from  six  points  of  the  5 1 .6  and  57  degree 
inclination  launches  generated  by  the  Gravtum  program.  From  each  of  the  six  points,  constant  pitch  and 
roll  angles  are  cycled  through  every  combination  between  30  and  90  degrees  at  5  degree  increments  for 
pitch  and  0  and  90  degrees  for  roll  and  numerically  integrated  up  to  the  nominal  landing  altitude  of  the 
orbiter.  At  this  point  the  longitude  and  latitude  is  plotted  against  the  map  of  the  Eastern  United  States 
creating  landing  zones  and  the  average  velocity  of  each  landing  zone  is  compared  to  the  nominal  Shuttle 
landing  speed.  Airport  locations  are  then  added  to  compare  intercept  points  for  possible  abort  landings. 

Results  drawn  from  the  data  determine  whether  the  constant  control  descent  proves  the  feasibility 
of  an  east  coast  abort  or  if  further  control  methods  are  necessary. 


3.2  Gravity-Turn  Trajectory  Launch 

Numerically  integrating  the  equations  described  in  2.2.3  and  plotting  the  latitude  and  longitude  at 
every  2  seconds  throughout  the  launch  gives  a  picture  of  the  launch  trajectory.  When  plotted  with  a 
computer  generated  map  of  the  Eastern  United  States  it  demonstrates  just  how  close  the  ascent  brings  the 
orbiter  to  the  east  coast.  Figures  3-1  and  3-2  show  the  launch  of  a  51 .6  and  57  degree  inclination  launch  out 
of  Cape  Canaveral  from  0  to  480  seconds. 
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Figure  3-1.  51,6  Degree  Inclination  Launch  from  Cape  Canaveral 


Figure  3-2.  57  Degree  Inclination  Launch  from  Cape  Canaveral 
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After  the  launch  was  plotted  six  points  throughout  the  launch  were  chosen  whose  state  vector  would  be 
used  as  the  initial  state  for  the  Abort  program.  Since  nothing  can  be  done  until  SRB  separation  the  first 
abort  point  was  selected  just  afterwards  at  120  seconds.  The  remaining  five  points  were  picked  at  60 
second  intervals  of  180, 240,  300, 360,  and  420  seconds.  Figures  3-3  and  3-4  show  these  six  abort  points 
mapped  to  the  51.6  and  57  degree  inclination  launches  shown  in  Figures  3-1  and  3-2. 
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3.3  Abort  Landing 

The  Abortl  program  uses  the  six  abort  points  as  initial  states  to  numerically  integrate  equations 
2.22  through  2.27.  Pitch  and  roll  angles  are  held  constant  throughout  the  descent,  but  cycled  through  all 
possible  combinations  ranging  from  30-90  degrees  of  pitch  and  0-90  degrees  of  roll  for  each  abort  point. 
The  ranges  in  pitch  exceed  the  angle  parameters  established  by  the  equations  in  the  aerodynamics  model 
discussed  in  2.3.4  but  provide  all  possibilities  within  the  Shuttle’s  flight  envelope.  This  range  increase 
simply  assumes  that  the  aerodynamics  equations  apply  for  increased  angles  of  45  degrees  in  the  hypersonic 
continuum  and  30  degrees  in  free  molecular  flow.  The  equations  are  integrated  using  all  these  angular 
combinations  until  an  altitude  of  25,000  meters  is  reached,  which  is  the  nominal  landing  altitude  of  the 
Orbiter.  The  final  latitude  and  longitude  of  each  combination  is  then  plotted  on  the  computer  generated 
map  of  the  Eastern  United  States  to  form  six  different  zones.  These  landing  zones  are  shown  in  figures  3-5 
and  3-6  for  both  a  5 1 .6  and  57  degree  inclination  launch. 
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Figure  3-5.  Abort  Landing  Zones  for  a  51 .6  Degree  Inclination  Launch 


Figure  3-6.  Abort  Landing  Zones  for  a  57  Degree  Inclination  Launch 
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These  abort  zones  represent  the  possible  locations  of  the  Obiter  at  25,000  meters  but  don’t  display 
the  speed.  For  the  altitude  of  25,000  meters,  the  nominal  approach  speed  of  the  Orbiter  is  Mach  2.5  (6: 1 ). 
The  speeds  of  each  combination  were  averaged  to  give  an  average  velocity  for  each  landing  zone.  Figures 
3-7  and  3-8  show  the  average  speed  of  each  zone  compared  to  the  nominal  approach  speed. 


Average  Velocity  of  Abort  Landing  Zones  for  a  51.6  Degree 
Inclination  Launch 
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Figure  3-7.  Average  Velocities  of  Abort  Landing  Zones  for  a  51.6  Degree  Inclination  Launch 


Average  Velocities  of  Abort  Landing  Zbnesfora  57  Degree 
Inclination  Launch 
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Figure  3-8.  Average  Velocities  of  Abort  Landing  Zones  for  a  57  Degree  Inclination  Launch 
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3.4  East  Coast  Airports 

Table  3-1  lists  thirty  seven  airports  close  to  the  launch  path  of  the  shuttle  citing  the  location  in 
longitude  and  latitude  and  giving  runway  lengths  in  feet . 


Table  3-L 


Airport 

Latitude 

Longitude 

Runway  Length  (ft) 

Andrews  AFB 

38.82 

-76.87 

9755 

Atlantic  City  Inti 

39.45 

-74.58 

10000 

Baltimore/Martin  State 

39.33 

-76.42 

6996 

Bangor  Inti 

44.8 

-68.83 

11439 

Beaufort  MCAS 

32.48 

-80.72 

12202 

Brunswick  NAS 

43.9 

-69.93 

8000 

Burlington  Inti 

44.47 

-73.15 

7821 

Cecil  Fid  NAS 

30.22 

-81.88 

12503 

Charleston  AFB 

32.9 

-80.03 

9001 

Cherry  Point 

34.9 

-76.88 

8980 

Dover  AFB 

39.13 

-75.47 

12902 

Gander 

48.93 

-54.57 

10500 

Goose  Bay 

53.02 

-60.42 

11000 

GrifRssAFLD 

43.23 

-75.4 

11820 

Langley  AFB 

37.08 

-76.35 

10000 

McEntire  ANGS 

33.92 

-80.8 

9001 

McGuire  AFB 

40.02 

-74.58 

10001 

Moncton 

46.17 

-64.57 

8000 

Narsarsuaq 

61.17 

-45.42 

6000 

Niagara  Falls  Inti 

43.1 

-78.95 

9125 

Norfolk  NAS 

36.93 

-76.28 

8369 

Oceana  NAS 

36.82 

-76.03 

11997 

OtisANGB 

41.65 

-70.52 

9500 

Patuxent  River  NAS 

38.28 

-76.42 

11800 

Pope  AFB 

35.17 

-79.02 

7501 

Richmond  Inti 

37.5 

-77.32 

8999 

Savannah  Inti 

32.13 

-81.2 

9351 

Seymour  Johnson  AFB 

35.33 

-77.97 

11758 

Shaw  AFB 

35.97 

-80.47 

10010 

Sondrestrom  AB 

67.02 

-50.72 

9200 

St.  Johns 

47.62 

-52.75 

8500 

Sydney 

46.17 

-60 

7100 

Syracuse/Hancock  Inti 

43.12 

-76.12 

9003 

West  Field/Bames  Muni 

42.17 

-72.72 

9000 

Westover  ARB 

42.2 

-72.53 

11600 

Willow  Grove  NAS 

40.2 

-75.15 

8002 

41.95 

-72.67 

9502 
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In  figures  3-9  and  3-10  these  airport  locations  were  plotted  onto  figures  3-5  and  3-6  to  determine 
which  would  be  the  most  likely  abort  landing  sites.  Primary  landing  sites  will  be  those  airports  within  the 
landing  zone  that  have  the  longest  runway  while  secondary  sites  will  be  within  or  close  to  the  zone.  It 
would  be  preferable  to  find  runways  of  10,000  ft.  or  more  that  could  safely  land  the  Orbiter  using  a  drogue 
chute  but  if  none  are  available  a  short  runway  is  much  better  than  ditching  into  the  ocean.  If  no  airports 
fall  within  the  zone  the  closest  airport  will  be  chosen  as  a  remote  possibility.  These  remote  sites  would 
mean  that  the  Orbiter  would  have  to  travel  outside  it’s  nominal  landing  approach  and  should  be  avoided  if 
possible. 
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Figure  3-10,  East  Coast  Airports  Compared  to  57  Degree  Inclination  Launch  Abort  Landing  Zones 
3.5  Results 

Tables  3-2and  3-3  list  possible  abort  landing  sites  for  both  51 .6  and  57  degree  inclination  launches 
using  the  guides  established  in  the  previous  section. 


Table  3-2.  Abort  Landing  Sites  for  51.6  Degree  Inclination  Launch 


Landing 

No.  of  Airports 

Primary 

Secondary 

Zone 

within  Zone 

Site 

Site 

1 

None 

**  Charleston  AFB 

None 

2 

None 

**  Cherry  Point 

None 

3 

None 

**  Oceana  NAS 

None 

4 

None 

**Ottis 

None 

5 

1 

Sydney 

Moncton 

6 

2 

Gander 

St.  Johns 
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Table  3-3.  Abort  Landing  Sites  for  57  Degree  Inclination  Launch 


Landing 

No.  of  Airports 

Primary 

Secondary 

Zone 

within  Zone 

Site 

Site 

1 

None 

**  Charleston  AFB 

None 

2 

None 

**  Cherry  Point 

None 

3 

None 

**  Oceana  NAS 

None 

4 

None 

**Ottis 

None 

5 

3 

Bangor 

Moncton 

6 

4 

Goose  Bay 

Gander 

From  the  tables,  it  is  plain  to  see  that  for  both  5 1 .6  and  57  degree  inclination  launches  that  the  only  airports 
that  lie  in  the  landing  zones  occur  at  360  and  420  seconds  into  the  launch.  At  these  later  times,  the  Orbiter 
is  high  enough  that  it  could  do  an  AOA  or  ATO  and  wouldn’t  need  to  abort  to  the  east  coast.  The  earlier 
abort  points  occurring  at  120  and  180  seconds  could  abort  using  a  RTLS  so  it  isn’t  imperative  to  try  a 
landing  at  Charleston  or  Cherry  Point.  The  midrange  landing  zones  however,  have  no  alternative  but  to 
attempt  the  abort  to  Oceana  and  Ottis,  which  are  far  outside  the  predicted  landing  zones  established  by  the 
Orbiter’s  nominal  landing  approach. 

Furthermore,  figures  3-7  and  3-8  indicate  that  the  velocity  of  the  orbiter  in  each  zone  is  slower 
than  the  nominal  mach  2.5  for  the  25,000  meter  altitude.  This  may  cause  the  Orbiter  to  change  it’s  glide 
slope 

in  order  to  make  up  for  this  reduction  in  energy.  If  the  speed  is  too  slow  it  may  not  make  the  runway  at  all 
and  park  the  Orbiter  in  someone’s  back  yard.  In  either  case  the  landing  zones  do  not  represent  a  normal 
shuttle  approach  for  landing. 

The  only  conclusion  that  can  be  draw  is  that  a  constant  control  descent  does  not  answer  the 
question  of  east  coast  abort  feasibility  but  only  provides  a  guide  for  fiirther  control  methods.  The  landing 
zones  2,  3,  and  4  for  both  inclination  launches  come  close  to  some  airports  whose  runways  could 
accommodate  the  Shuttle .  These  airports  now  become  selected  landing  sites  which  the  second  control 
descent  program,  Abort2,  will  attempt  to  reach. 
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4.  Data  Analysis  and  Results  for  Controlled  Descent  to  Landing  Site 


4. 1  Introduction 

The  controlled  descent  method  picks  primaiy  airports  based  on  location  and  runway  size  and 
guides  the  orbiter  to  them.  Pitch  and  roll  only  stay  constant  until  the  correct  heading  angle  to  intercept  the 
primary  airport  is  reached,  then  the  roll  angle  becomes  zero  allowing  the  orbiter  to  glide  straight  toward  the 
target.  This  is  done  for  every  point  throughout  the  trajectory  for  both  5 1 .6  and  57  degree  inclination 
launches  giving  exact  times  when  primaiy  aiiports  can  be  reached.  Average  velocity  and  altitude  are  given 
for  all  abort  attempts  within  a  zone  around  the  target  airport  which  is  1  degree  of  longitude  by  1  degree  of 
latitude  with  the  runway  in  the  center.  These  averages  are  then  compared  with  the  nominal  landing  criteria 
of  a  speed  of  mach  land  an  estimated  altitude  of 6098  meters  (6:1). 

4.2  Controlled  Descent  Using  Heading  Angle 

The  Gravtum  program  once  again  simulates  the  launch  of  the  shuttle  at  both  5 1 .6  and  57  degrees 
but  the  abort  program  is  slightly  modified.  Rather  than  cycling  through  all  possible  angles  to  create  landing 
zones  for  six  discrete  points,  an  airport  is  chosen  as  a  landing  site  and  every  point  in  the  laimch  from  120  to 
480  seconds  at  2  second  intervals  is  aborted  to  that  site.  This  is  done  by  choosing  both  pitch  and  roll  as  45 
degrees  which  are  reasonable  values  that  are  within  the  flight  envelope  established  in  the  Shuttle’s 
aerodynamic  flight  model.  These  conditions  remain  constant  until  the  proper  heading  angle  to  intercept  the 
target  airport  is  achieved  at  wiiich  point  the  roll  and  pitch  angle  become  0  and  30  degrees  respectively. 

Spherical  trigonometry  is  used  to  calculate  the  proper  heading  angle  to  the  prescribed  airport 
knowing  the  latitude  and  longitude  of  the  target  as  well  as  the  latitude  and  longitude  of  the  Orbiter  at  every 
point.  The  following  two  equations  were  derived  from  spherical  trigonometry,  because  you  must  first  solve 
one  to  obtain  the  other. 
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where 

/  =  Radial  Distance  to  Target  (radians) 
xjfAP  ~  Heading  Angle  to  Airport  (radians) 
Sap  =  Latitude  of  Airport  (radians) 

SoR  =  Latitude  of  Orbiter  (radians) 

=  Longitude  of  Airport  (radians) 

XoR  =  Longitude  of  OftzVer(radians) 


These  two  equations  were  placed  in  the  Dynam  subroutine  of  the  Abort2  program  so  that  each  abort  descent 
would  roll  to  0  degrees  at  the  proper  heading  angle. 

4.3  Abort  Landing  Sites 

From  the  abort  landing  zones  established  in  the  constant  control  method  certain  airports  were 
chosen  as  ‘primary’  and  ‘secondary’  landing  sites.  Due  to  the  lack  of  airports  within  the  early  abort  regions 
several  other  airports  were  chosen  as  ‘possible’  landing  sites  if  alternate  control  methods  were  used.  For 
both  the  51.6  and  57  degree  inclination  laimches  the  same  landing  sites  for  the  early  aborts  were  chosen 
because  neither  laimch  produced  landing  zones  which  intercepted  any  airports.  Since  both  launches  have 
identical  landing  sites  for  the  early  abort  points,  the  57  degree  launch  is  investigated  first  because  it’s 
trajectory  brings  it  closer  to  the  east  coast  and  if  the  airports  cannot  be  reached  there  is  no  point  in  doing  the 
abort  from  the  5 1 .6  degree  laimch. 

The  first  three  landing  sites  chosen  were  Charleston  AFB,  Cherry  Point,  and  Oceana  Naval  Air 
Station  (NAS).  These  sites  represent  the  first  three  minutes  of  abort  opportunity  starting  after  SRB 
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separation  at  120  seconds  and  ending  at  300  seconds.  Figures  4-1, 4-2,  and  4-3  show  the  abort  descents 
from  the  launch  trajectory  down  to  6098  meters. 


Figure  4-1 .  Abort  Landing  Attempt  at  Charleston  AFB 


Figure  4-2.  Abort  Landing  Attempt  at  Cherry  Point 
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Figure  4-3.  Abort  Landing  Attempt  at  Oceana  NAS 

These  figures  clearly  illustrate  that  none  of  these  three  airports  can  be  approached  under  nominal  conditions 
using  a  non  powered  control  descent  from  a  57  degree  inclination  launch  with  pitch  and  roll  conditions  as 
mentioned  in  section  4.2. 
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The  end  of  each  abort  represents  the  estimated  nominal  altitude  of  6098  meters  so  it  may  be  possible  to 
reach  these  airports  but  the  altitude  would  be  too  low  to  allow  sufficient  room  to  align  with  the  runway  and 
land  safely.  These  results  make  running  the  abort  descents  for  a  5 1 .6  degree  inclination  unnecessary  due  to 
the  longer  distances  that  would  need  to  be  covered. 

Ottis  Air  National  Guard  Base  (ANGB)  was  the  next  landing  site  chosen  for  both  57  and  5 1 .6 
degree  inclination  launches.  This  site  covers  the  abort  opportunities  from  300  to  360  seconds  into  the 
launch.  Figures  4-4  shows  that  there  are  intercepts  for  the  57  degree  laxmch  but  the  window  of  opportunity 
is  only  8  seconds  at  330  seconds  into  the  launch.  Figure  4-5  shows  that  no  intercepts  are  possible  for  the 
51.6  degree  inclination.  These  outcomes  leave  the  51.6  degree  launch  with  no  landing  opportunities  and 
the  57  degree  launch  with  a  very  small  one. 

The  next  landing  sites  chosen  were  Sydney  for  the  51.6  degree  launch  and  Moncton  for  the  57 
degree  inclination  launch.  The  reason  for  choosing  Moncton  over  Bangor  was  the  result  of  reducing  the 
flight  envelope  from  90  degree  pitch  to  45  degree  pitch.  This  reduction  doesn’t  allow  the  Orbiter  to  slow 
down  fast  enough  to  reach  Bangor  so  the  secondary  site  of  Moncton  was  made  the  primary.  Figures  4-6  and 
4-7  illustrate  these  abort  landings. 
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Figure  4-6.  Abort  Landing  at  Moncton  for  a  57  Degree  Inclination  Launch. 


Figure  4-7.  Abort  Landing  at  Sydney  for  a  5 1 .6  Degree  Inclination  Launch 
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Figure  4-6  proves  that  Moncton  is  within  reach  for  many  descents  starting  at  362  seconds  into  the  launch  up 
until  396  seconds.  Figure  4-7  on  the  other  hand  shows  that  an  abort  to  Sydney  is  possible  but  the  window 
of  opportunity  is  much  smaller  at  12  seconds  starting  at  390  seconds  into  the  5 1.6  degree  inclination  launch. 

The  final  abort  airports  Gander,  St.  Johns  and  Goose  Bay  that  represent  420  to  480  seconds  into 
the  launch  were  too  close  in  the  landing  zone  to  be  reachable  using  a  45  degree  pitch.  The  remaining 
airports  on  the  rim  of  the  landing  zones  like  Narsarsuaq  and  Sondrestrom  were  out  of  reach.  Different 
configurations  of  pitch  and  roll  which  are  out  of  the  flight  envelope  may  allow  the  Orbiter  to  reach  these 
sites  but  velocity  and  altitude  will  be  vital  in  determining  whether  the  Orbiter  can  bum  off  enough  energy  to 
become  slow  enough  and  low  enough  to  land. 

4, 4  Altitude  and  Velocity  Profiles 

The  velocities  and  altitudes  of  all  the  descent  trajectories  within  a  1  by  1  degree  zone  around  the 
landing  sites  were  averaged.  Figures  4-8  and  4-9  compare  these  averages  to  the  nominal  approach  velocity 
and  altitude  of  mach  1  at  6098  meters.  The  speed  of  soimd  at  6098  meters  was  computed  to  be  347  m/s 
using  the  subroutine  ATM  and  used  as  the  nominal  speed  of  approach. 


Average  Velocities  Compared 
to  Nominal  Approach  Velocities 


Figure  4-8.  Average  Velocities  for  All  Airport  Intercepts 
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Average  Altitudes  Compared 
to  Nominal  Approach  Altitude 
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Figure  4-9.  Average  Altitudes  for  All  Airport  Intercepts 


From  these  charts  Sydney  and  Moncton  are  too  high  and  too  fast  where  as  Ottis  is  too  high  and  not  fast 
enough.  These  results  don’t  mean  the  Orbiter  can’t  land  at  these  sites  but  what  they  do  suggest  is  that 
energy  management  techniques  such  as  "S’  turns  and  high  pitch  angles  may  need  to  be  employed  to  reduce 
speed  and  in  turn  altitude.  To  increase  speed  and  reduce  altitude  a  lower  pitch  angle  might  be  employed  . 

4, 5  Results  and  Recommendations 

The  method  of  controlled  descent  using  heading  angle  for  a  5 1 .6  and  57  degrees  inclination  launch 
yield  abort  procedures  outlined  in  table  4-1  and  4-2. 


Table  4-1 .  Abort  Procedures  for  a  5 1 .6  Degree  Inclination  Launch 


Time  into  Launch  (seconds) 

Abort  Mode 

120-280 

RTLS 

281-389 

Contingency 

390-402 

ECA  to  Sydney 

403  -  up 

Contingency,  AO  A,  or  ATO 
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Table  4-2>  Abort  Procedures  for  a  57  Degree  Inclination  Launch 


Time  into  Launch  (seconds) 

Abort  Mode 

120-280 

RTLS 

281-329 

Contingency 

330-338 

ECA  to  Ottis  ANGB 

339-361 

Contingency 

362-396 

ECA  to  Moncton 

397-up 

Contingency,  AOA,  or  ATO 

As  shown  in  both  tables  the  Return  To  Launch  Site  (RTLS)  abort  mode  must  still  be  used  up  to  280 
seconds  into  the  launch.  This  is  due  to  the  fact  that  the  Orbiter  is  not  high  or  fast  enough  to  reach  any  of  the 
airports  in  the  Southeast.  For  the  51 .6  degree  case  the  next  110  seconds  would  be  a  contingency  or  ditch 
type  of  abort  as  would  the  57  degree  case  but  only  lasting  50  seconds.  The  51.6  degree  case  can  then  East 
Coast  Abort  (EC  A)  to  Sydney  for  12  seconds  after  which  depending  on  engine  capability  the  Orbiter  could 
ditch,  Abort  Once  Around  (AO A),  or  Abort  To  Orbit  (ATO).  The  57  degree  case  has  two  more 
opportumties  to  EC  A  to  Ottis  and  Mwicton  at  330  and  362  respectively  before  having  to  either  ditch,  AOA, 
or  ATO. 

This  investigation  has  shown  that  ECA’s  are  possible  but  very  limited.  The  early  abort  cases  are 
still  left  up  to  the  RTLS  and  contingency  modes  which  are  either  extremely  experimental  or  result  in  the 
loss  of  the  Orbiter.  Further  study  into  ECAs  using  optimal  controls  may  yield  larger  windows  of 
opportumty  for  the  abort  landing  sites  mentioned  as  well  as  reduce  speed  and  velocity  to  ensure  safe 
landings.  As  for  the  early  abort  cases  and  the  51 .6  degree  launches,  a  powered  case  should  also  be 
investigated  to  maneuver  the  orbiter  closer  to  the  coast  line  before  releasing  the  external  tank .  This  may 
allow  the  51.6  degree  launches  to  reach  the  landing  sites  that  the  non-powered  case  could  not  as  well  as 
give  both  57  and  51.6  launches  the  flexibility  of  choosing  landing  sites  that  may  have  longer  runways  and 
better  facilities  but  were  too  far  inland  to  reach. 
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In  conclusion  this  thesis  does  indeed  prove  that  East  Coast  Aborts  are  possible  but  more 
importantly  proves  that  further  research  should  be  done  to  improve  landing  opportunities  as  well  as  increase 
the  number  of  reachable  airports  in  which  to  land.  This  study  must  also  be  followed  up  to  investigate  the 
aerodynamic  moments  the  Shuttle  may  experience  as  a  result  of  hypersonic  maneuvers  which  were  not 
looked  at  in  this  point  mass  model.  In  any  case  future  research  is  the  major  recommendation  of  this  paper 
because  with  the  greater  number  of  launches  the  shuttle  makes  the  larger  the  chances  that  something  could 
go  wrong  and  to  be  able  to  save  both  crew  and  Orbiter  would  be  beneficial.  Not  only  does  this  benefit 
NASA  but,  if  we  can  avoid  a  major  catastrophe  and  prevent  further  budget  cuts  or  an  all  out  cancellation  of 
the  space  program,  then  it  would  also  benefit  everyone  who  has  ever  looked  up  to  the  stars  and  dreamed. 
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Appendix  A  Gravity-Turn  Trajectory  Launch  Program 


Contained  in  this  section  is  all  the  source  code  used  in  modeling  the  launch  of  the  Space  Shuttle. 
The  only  subroutines  that  are  omitted  are  the  two  Haming  numerical  integrators  use  to  integrate  Rhs2  and 
Rhs3.  They  are  left  out  because  the  code  is  the  same  except  for  which  Rhs  subroutine  it  calls  to  integrate. 


Program  Gravtum 

implicit  none 
c 

c  dynamics  model  for  gravity-turn  trajectory  fi’om  Wiesel  page  208  equations  67  -  78. 
c 

c  Input:  Initial  State  ( Down  Range  Distance,  Altitude,  Velocity,  Flight  Path  Angle, Mass) 
c  Output:  State  vector  at  any  time  T  ( X,  H,  V,  Gamma) 
c  (Printed  to  screen  and  selected  files) 

c 

common  /ham/  t,x(230,4),f(230,4),err(230),n,h 
Real*  8 1  ,Long_Launch,Lat_Launch,Inclin,R 
Real*8  x  ,  f ,  err ,  h ,  toler,  Long  ,  Lat,  Psi 
integer  nxt,  n 
toler  =  l.Od-9 


c  File  used  in  controlled  Abort 

open  (unit=  2, file  =  'c:\afit\grav.out’,  status  -  unknown’) 

c  File  used  in  mapping  program  for  ascent  trajectory 

open  (imit=  3,file  =  'c:\afit\grav2. out’,status  =’unknown') 

c  Files  used  in  constant  control  abort  and  in  the  mapping  program 

open  (unit=4,  file  =  ’c:\afit\abortl  .out', status  -untoown') 
open  (unit==8,  file  =  'c:\afit\abort2.out',status  -  unknown') 
open  (imit=9,  file  =  'c:\afit\abort3  .out’,status  -  unknown') 
open  (unit=10,  file  =  'c:\afit\abort4.out', status  -  unknown') 
open  (unit=l  1,  file  =  'c:\afit\abort5.out’,status  -unknown') 
open  (unit=12,  file  =  'c:\afit\abort6.out', status  ='unknown') 
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initialize(tiine, state)  and  constants  (entered  as  SI  units) 


x(l,l)  =  0 

x(2,l)=100 

x(3,l)  =  32.1 

x(4,l)  =  1.56209405 

x(5,l)  =  2041200 

Long_Launch  =  -80.6 

LatLaunch  =  28.4583 

Inclin  =  51.6 

R  =  6378.135d3 

n=5 

t=0.d0 

h=2.0d0 

nxt=0 

call  hamingl(nxt) 
if(nxt.eq.0)then 
write(*,*)  ’dead’ 
stop 
endif 

First  stage  on  launch  profile 
do  while(t.le.  120.and.x(4,nxt).gt.toler) 

call  hamingl(nxt) 

Calculations  for  Long,  Lat,  and  Psi 

Lat  =  -dacosd(dcosd(90-Lat_Launch)*dcos(x(l,nxt)/R) 

1  +  dsind(90'Lat_Launch)*dsin(x(l  ,nxt)/R) 

2  *dcosd(90-Inclin))+90 

Long  =  dacosd((dcos(x(l,nxt)/R)-dcosd(90-Lat) 

1  *dcosd(90“Lat_Launch))/(dsind(90-Lat) 

2  *dsind(90“Lat_Launch)))4-  Long_Launch 

Psi  =1.5708  -  (dasin(dcosd(Inclin)/dcosd(Lat))) 
write(3,*)  Long, Lat 

write(*,*)'state  =',  x(l,nxt),x(2,nxt),x(3,nxt),  x(4,nxt), 

1  x(5,nxt),t 

if(t.eq,120)then 

write(4,*)  x(2,nxt)+R,Long,  Lat,x(3,nxt),x(4,nxt),Psi 
For  state  at  abort  point  1 

write(2,*)  x(2,nxt)+R,Long,  Lat,x(3,nxt),x(4,nxt),Psi 
endif 
end  do 

Initializing  for  second  stage 
x(l,l)=  x(l,nxt) 

x(2,l)=  x(2,nxt) 
x(3,l)=  x(3,nxt) 
x(4,l)=  x(4,nxt) 
x(5,l)=  x(5,nxt)- 164000 


call  haining2(nxt) 
if(nxt.eq.O)then 

write(*,*)  ’dead2' 

stop 

else 

write(*,*)  ’alive’ 
endif 

second  stage  of  launch  profile 
do  while(t.lt.480) 

call  haming2(nxt) 

Calculations  for  Long,  Lat,  and  Psi 

Lat  =  -dacosd(dcosd(90-Lat_Launch)*dcos(x(l,nxt)/R) 

1  +  dsind(90-Lat_Launch)*dsin(x(  1  ,nxt)/R) 

2  *dcosd(90"Inclin))+90 

Long  =  dacosd((dcos(x(l,nxt)/R)-dcosd(90-Lat) 

1  *dcosd(90-Lat_Launch))/(dsind(90-Lat) 

2  *dsind(90-Lat_Launch)))+  Long_Launch 

Psi  =  1.5708  “  (dasin(dcosd(Inclin)/dcosd(Lat))) 
vvrite(3,*)  Long, Lat 

write(*,*)'state  =’ ,  x(l,nxt),x(2,nxt),x(3,nxtXx(4,nxt), 

1  x(5,nxt),t 

write(2,*)  x(2,nxt)+R,Long,  Lat,x(3,nxt),x(4,nxt),Psi 

For  abort  point  2 
if(t.eq.l80)then 

write(8,*)  x(2,nxt)+R,Long,  Lat,x(3,nxt),x(4,nxt),Psi 
endif 

For  abort  point  3 
if(t.eq.240)then 

write(9,*)  x(2,nxt)+R,Long,  Lat,x(3,nxt),x(4,nxt),Psi 
endif 

For  abort  point  4 
if(t.eq.300)then 

write(10,*)  x(2,nxt)+R,Long,  Lat,x(3,nxt),x(4,nxt),Psi 
endif 

For  abort  point  5 
if(t.eq.360)then 

write(l  1,*)  x(2,nxt)+R,Long,  Lat,x(3,nxt),x(4,nxt),Psi 
endif 

For  abort  point  6 
if(t.eq.420)then 

write(12,*)  x(2,nxt)-4-R,Long,  Lat,x(3,nxt),x(4,nxt),Psi 
endif 


c  initialize  for  final  stage 

x(l,l)=  x(l,nxt) 
x(2,l)=  x(2,nxt) 
x(3,l)=  x(3,nxt) 
x(4,l)=  x(4,nxt) 
x(5,l)=  x(5,nxt)-32000 
nxt=0 

call  haming3(nxt) 
if(nxt.eq.O)then 
write(*,*)  ’dead3' 
stop 
endif 

c  Final  stage 

do  while(t.ge.480) 

call  haming3(nxt) 

c  Calculations  for  Long,  Lat,  and  Psi 

Lat  =  -dacosd(dcosd(90“Lat_Laiinch)*dcos(x(l,nxt)/R) 

1  +  dsind(90-Lat_Launch)*dsin(x(l  ,nxt)/R) 

2  *  dcosd(90-Inclin))+90 

c 

Long  =  dacosd((dcos(x(l,nxt)/R)-dcosd(90-Lat) 

1  *dcosd(90-Lat_Launch))/(dsind(90-Lat) 

2  *dsind(90-Lat_Launch)))-f  LongLaunch 
c 

Psi  =  1.5708  -  (dasin(dcosd(Inclin)/dcosd(Lat))) 
c 

write(3,*)  Long, Lat 

write(2,*)  x(2,nxt)+R,Long,  Lat,x(3,nxt),x(4,nxt),Psi 
write(*,*)'state  x(l,nxt),x(2,nxt),x(3,nxt),x(4,nxt), 

1  x(5,nxt),t 

stop 
end  do 

write(*,*)  Tinished  Successfully  with' 
write(*,*) '  Gamma  = x(4,nxt),  'at  time',t  ,nxt 
close(2) 
close(3) 
close(4) 
close(8) 
close(9) 
close(lO) 
close(ll) 
close(12) 
stop 

end 
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Subroutine  Rhsl(nxt) 


c 

c  Equations  of  Motion  for  gravity-turn  trajectory  from  Wiesel  page  208  equations  67  -  78. 

c  covering  from  t=0s  to  t=  1 20s 

c 
c 

c  RHS  needs  access  to  the  common  /ham/: 

common  /ham/  t,x(230,4),  f(230,4),err(230),n,h,mode 
real*8  t,x,f,err,h,mu,R,m_dotl  ,Tr,P0,PALT,TALT,DALT,dDdr 
real*8  CDAl 

c  List  of  Constants 

R=6378.135d3 
mu  =3.98601dl4 
mdotl  =9906 
Tr=30300d3 
P0=  101325.0d0 
CDAl  =  101.12 


c  Must  include  the  gravity  term 

c 

g  =  mu/((x(2,nxt)+  R)*(x(2,nxt)+  R)) 
c 

c  Equations  of  Motion  from  Wiesel 

c 

c  f(l)  -  d(x)/dt  -  Down  Range  Distance 
c  f(2)-d(H)/dt- Altitude 
c  f(3)  -  d(V)/dt  -  Velocity 
c  f(4)  -  d(gamma)/dt  -  Flight  Path  Angle 
c  f(5)  -  d(mass)/dt  -  Mass 
c 

c  Wiesel  Equations  67-70  (including  Drag) 
c 

f(l,nxt)=  x(3,nxt)*dcos(x(4,nxt)) 

f(2,nxt)=  x(3,nxt)*dsin(x(4,nxt)) 

c  Calling  subroutine  ATM  to  get  Density  at  Altitude 

callATM(x(2,nxt),P0,PALT,TALT,DALT,dDdr) 

f(3,nxt)=  (Tr/x(5,nxt))  -  (g*dsin(x(4,nxt))) 

1  +  ((f(l,nxt)*f(l,nxt))/(x(2,nxt)+R))*dsin(x(4,nxt)) 

2  -  ((CDAl*.5*DALT*x(3,nxt)*x(3,nxt))/x(5,nxt)) 

f(4,nxt)=-(g*dcos(x(4,nxt))/(x(3,nxt))) 

1  +  (f(l,nxt)*f(l^t)*dcos(x(4,nxt))/(x(3,nxt)*(x(2,nxt)+R))) 

f(5,nxt)=  -mdotl 

return 

end 
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Subroutine  Rhs2(nxt) 

c 

c  Equations  of  Motion  for  gravity-turn  trajectory  from  Wiesel  page  208  equations  67  -  78. 

c  For  second  stage  from  120s  to  480s 
c 
c 

c  RHS  needs  access  to  the  common  /ham/; 

common  /ham/  t,x(230,4),  f(230,4),err(230),n,h,mode 
real*8t,x,f,err,h,mu,R,m_dot2,Tr2,P0,PALTJALT,DALT,dDdr 


R=6378.135d3 
mu  =3.98601dl4 
m_dot2  =1496.0d0 
c  frill  Thrust 

c  Tr2  =6300d3 

c  partial  Thrust 

Tr2  =  6000d3 

CDA2  =  79.5 
P0=101325.0d0 


c  Must  include  the  gravity  term 

c 

g  =  mu/((x(2,nxt)+  R)*(x(2,nxt)+  R)) 
c 

c  Equations  of  Motion  from  Wiesel 

c 

c  f(  1 )  -  d(x)/dt  -  Down  Range  Distance 

c  f(2)-d(H)/dt- Altitude 
c  f(3)-d(V)/dt- Velocity 
c  f(4)  -  d(gamma)/dt  -  Flight  Path  Angle 
c  f(5)  -  d(mass)/dt  -  Mass 
c 

c  Wiesel  Equations  67-70  (Including  Drag) 
c 

f(l,nxt)=  x(3,nxt)*dcos(x(4,nxt)) 
f(2,nxt)=  x(3,nxt)*dsin(x(4,nxt)) 
c  Calling  subroutine  ATM  to  get  Density  at  Altitude 

callATM(x(2,nxt),P0,PALT,TALT,DALT,dDdr) 

f(3,nxt)=  (Tr2/x(5,nxt))  -  (g*dsin(x(4,nxt))) 

1  +  ((f(l,nxt)’'‘f(l,nxt))/(x(2,nxt)+R))*dsin(x(4,nxt)) 

2  -  ((CDA2*.5*DALT*x(3,nxt)*x(3,nxt))/x(5,nxt)) 
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c  To  control  shape  of  ascent  by  positive  input  to  gamma 

if(t.eq.210)then 
x(4,nxt)  =  x(4,nxt)  +  .11 
endif 

if(t.eq.280)then 
x(4,nxt)  =  x(4,nxt)  +  .12 
endif 

if(t.eq.340)then 
x(4,nxt)  =  x(4,nxt)  +  .10 
endif 

if(t.eq.420)then 
x(4,nxt)  =  x(4,nxt)  +  .02 
endif 

f(4,nxt)=  -(g*dcos(x(4,nxt))/(x(3,nxt))) 

1  +  (f(l,nxt)*f(l,nxt)*dcos(x(4,nxt))/(x(3,nxt)*(x(2,nxt)+R))) 

f(5,nxt)=  -m_dot2 

return 

end 
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Subroutine  Rhs3(nxt) 
c 

c  Equations  of  Motion  for  gravity-tum  trajectory  from  Wiesel  page  208  equations  67  -  78. 
c  For  final  stage  for  over  480s 
c 
c 

c  RHS  needs  access  to  the  common  /ham/: 

common /ham/  t,x(230,4),  f(230,4),erT(230),n,h,mode 
real*8t,x,f,err,h,mu,R,m_dot3,Tr3,P0,PALT,TALT,DALT,dDdr 


R=6378.135d3 
mu=3.98601dl4 
m_dot3  =0 
Tr3=0 

P0=  101325.0d0 
CDA3  -  29.235 


c  Must  include  the  gravity  term 

c 

g  =  mu/((x(2,nxt)4-  R)*(x(2,nxt)+  R)) 
c 

c  Equations  of  Motion  from  Wiesel 

c 

c  f(l)  -  d(x)/dt  -  Down  Range  Distance 

c  f(2)-d(H)/dt~  Altitude 

c  f(3)-d(V)/dt- Velocity 
c  f(4)  -  d(gamma)/dt  -  Flight  Path  Angle 
c  f(5)  -  d(mass)/dt  -  Mass 
c 

c  Wiesel  Equations  67-70  (Including  Drag) 
c 

f(l,nxt)=  x(3,nxt)*dcos(x(4,nxt)) 
f(2,nxt)=  x(3,nxt)*dsin(x(4,nxt)) 
c  Calling  subroutine  ATM  to  get  Density  at  Altitude 

call  ATM(x(2,nxt),P0,PALT,TALT,DALT,dDdr) 
c 

f(3,nxt)=  (Tr3/x(55nxt))  -  (g*dsin(x(4,nxt))) 

1  +  ((f(l,nxt)*f(l,nxt))/(x(2,nxt)+R))*dsin(x(4,nxt)) 

2  -  ((CDA3*.5*DALT*x(3,nxt)*x(3,nxt))/x(5,nxt)) 

f(4,nxt)=  -(g*  dcos(x(4,nxt))/(x(3  ,nxt))) 

1  +  (f(  1  ,nxt)*f(l  ,nxt)*dcos(x(4,nxt))/(x(3,nxt)*(x(2,nxt)+R))) 

f(5,nxt)=  -m_dot3 

return 

end 
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Subroutine  Hamingl(nxt) 

c 

c  haming  is  an  ordinaiy  differential  equations  integrator 
c  it  is  a  fourth  order  predictor-corrector  algorithm 
c  which  means  that  it  carries  along  the  last  four 
c  values  of  the  state  vector,  and  extrapolates  these 
c  values  to  obtain  the  next  value  (the  prediction  part) 
c  and  then  corrects  the  extrapolated  value  to  find  a 
c  new  value  for  the  state  vector, 
c 

c  the  value  nxt  in  the  call  specifies  which  of  the  4  values 
c  of  the  state  vector  is  the  "next”  one. 
c  nxt  is  updated  by  haming  automatically,  and  is  zero  on 
c  the  first  call 
c 

c  the  user  supplies  an  external  routine  rhs(nxt)  which 
c  evaluates  the  equations  of  motion 
c 

common  /ham/  x,y(230,4),f(230,4),errest(230),n,h,mode 
double  precision  x,y,f,errest,h,hh,xo,tol,eps 
c 

c  all  of  the  good  stuff  is  in  this  common  block, 
c  X  is  the  independent  variable  ( time ) 
c  y(6,4)  is  the  state  vector-  4  copies  of  it,  with  nxt 
c  pointing  at  the  next  one 

c  f(6,4)  are  the  equations  of  motion,  again  four  copies 
c  a  call  to  rhs(nxt)  updates  an  entry  in  f 

c  errest  is  an  estimate  of  the  truncation  error  -  normally  not 
c  used 

c  n  is  the  number  of  equations  being  integrated  -  6  or  42  here 
c  h  is  the  time  step 

c  mode  is  0  for  just  EOM,  1  for  both  EOM  and  EOV 
c 

tol  =  O.OOOOOOOld+00 

c  switch  on  starting  algorithm  or  normal  propagation 
if(nxt)  190,10,200 
c 

c  this  is  hamings  starting  algorithm... .a  predictor  -  corrector 
c  needs  4  values  of  the  state  vector,  and  you  only  have  one-  the 
c  initial  conditions. 

c  haming  uses  a  Picard  iteration  (slow  and  painfull)  to  get  the 
c  other  three. 

c  if  it  fails,  nxt  will  still  be  zero  upon  exit,  otherwise 
c  nxt  will  be  1 ,  and  you  are  all  set  to  go 
c 

10xo  =  x 
hh  =  h/2.0d+00 
call  rhsl(l) 
do  401  =  2,4 
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x  =  x  +  hh 
do  20  i  =  l,n 

20  y(i,l)  =  y(i,l-l)  +  hh*f(i,I-l) 
call  rhsl(l) 
x  =  x  +  hh 
do30i=  l,n 

30y(U)  =  y(U.l)  +  h*f(i,l) 

40  call  rhsl(l) 
jsw  =  “20 
50  isw  =  1 
do  120  i  =  l,n 

hh  =  y(i,l)  +  h*(  9.0d+00*f(i,l)  +  19.0d+00*f(i,2) 

1  -  5.0d+00*f(i,3)  +  f(i,4) )  /  24.0d+00 

if(y(i,2)  .eq.  O.dO)  then 
eps  =  dabs(  hh  -  y(i,2) ) 
else 

eps  =  dabs(  (hh  -  y(i,2))/y(i,2) ) 
endif 

if(  eps  .It.  tol )  go  to  70 
isw  =  0 
70  y(i52)  =  hh 

hh  =  y(i,l)  +  h*(  f(i,l)  +  4.0d+00*f(i,2)  +  f(i,3))/3.0d+00 
if(y(i,3)  .eq.  O.dO)  then 
eps  =  dabs(  hh  -  y(i,3) ) 
else 

eps  =  dabs(  (hh-y(i,3))/y(i,3) ) 
endif 

if(  eps  .It.  tol )  go  to  90 
isw  =  0 
90y(i,3)  =  hh 

hh  =  y(i,l)  +  h*(  3.0d+00*f(i,l)  +  9.0d+00*f(i,2)  +  9.0d+00*f(i,3) 
1  +  3.0d+00*f(i,4) )  /  8.0d+00 

if(y(i,4)  .eq.  O.dO)  then 
eps  =  dabs(  hh  -  y(i,4)  ) 
else 

eps  =  dabs(  (hh-y(i,4))/y(i,4) ) 
endif 

if(  eps  .It.  tol )  go  to  1 1 0 
isw  =  0 

110y(i,4)  =  hh 
120  continue 
x  =  xo 

do  130  1  =  2,4 
x  =  x  +  h 
130  callrhsl(l) 
if(isw)  140,140,150 
140  jsw = jsw  +  1 
if(jsw)  50,280,280 
150x  =  xo 
isw=  1 
jsw=  1 
do  160  i  =  l,n 
160  eiTest(i)  =  0.0 
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nxt  =  1 
go  to  280 
190  jsw  =  2 
nxt  =  iabs(nxt) 

c  this  is  hamings  normal  propagation  loop  - 
c 

200  X  =  X  +  h 

npl  =  mod(nxt,4)  +  1 
go  to  (210,230),isw 
c  permute  the  index  nxt  modulo  4 
210  go  to  (270,270, 270, 220), nxt 
220  isw  =  2 

230  nm2  =  mod(npl  ,4)  +  1 
nml  =  mod(nm2,4)  +  1 
npo  =  mod(nml  ,4)  +  1 
c 

c  this  is  the  predictor  part 
c 

do  240  i  =  l,n 

f(i,nm2)  =  y(i,npl)  +  4.0d+00*h*(  2.0d+00*f(i,npo)  -  f(i,nml) 

1  +  2.0d+00*f(i,nm2) )  /  3.0d+00 

240  y(i,npl)  =  f(i,nm2)  -  0.925619835d0*errest{i) 
c 

c  now  the  corrector  -  fix  up  the  extrapolated  state 
c  based  on  the  better  value  of  the  equations  of  motion 
c 

call  rhsl(npl) 
do  250  i  =  l,n 

y(i,npl)  =  ( 9.0d+00*y(i,npo)  -  y(i,nm2)  +  3.0d+00*h*(  f(i,npl) 
1  +  2.0d+00*f(i,npo)  -  f(i,nml) ) )  /  8.0d+00 

errest(i)  =  f(i,nm2)  -  y(i,npl) 

250  y(i,npl)  =  y(i,npl)  +  0.0743801653d0  *  errest(i) 
go  to  (260,270) Jsw 
260  callrhsl(npl) 

270  nxt  =  npl 
280  return 
end 
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Subroutine  Atiii(aIt,pO,palt,talt,daIt,dddr) 

c 

c  earth  atmosphere  program,  regan  and  anandarskarian,  aiaa 
c  "dynamics  of  atmospheric  re-entiy",  appendix  a 
c 

c  input:  alt  altitude  in  meters 

c  po  ground  level  pressure,  n/sq  m 

c  output:  palt  pressure  at  altitude,  n/sq  m 
c  talt  temperature  at  altutude,  deg  c 

c  dalt  density  at  altitude,  kg/cu  m 

c  dddr  density  gradient,  kg/mM 

c 

double  precision  z(2 1  ),tm(2 1  ),lr(2 1  ),b,g0,r,d(2 1  ),p(2 1  ),pO,dO,rr 
double  precision  talt,palt,dalt,alt,m(21),dd(21) 
double  precision  el,e2,e3,e4,e5,re,galt 
double  precision  dddr,deldr,de2dr,de3dr 
c 

c  data  stmts  for  break  altitudes,  temperatures,  and  molecular  wts 
c 

c  altitudes 

data(z(i),i=l,21)/0.d3,  11.0191d3,  20.063 ld3,  32.1619d3, 

1  47.3501d3,  51.4125d3,  71.8020d3,  86.00d3, 

2  100.d3,  110.d3,  120.d3,  150.d3, 

3  160.d3,  170.d3,  190.d3,  230.d3, 

4  300.d3,  400.d3,  500.d3,  600.d3, 

5  700.d3  / 

c  molecular  temperature 

data(tm(i),i=l,21)/300.d0,  216.65d0,  216.65d0,  228.65d0, 

1  270.65d0,  270.65d0,  2]4.65d0,  ]86.946d0, 

2  210.65d0,  260.65d0,  360.65d0,  960.65d0, 

3  1110.60d0, 1210.65d0,  1350.65d0,  1550.65d0, 

4  1830.65d0, 2160.65d0,  2420.65d0, 2590.65d0, 

5  2700.0d0/ 
c  molecular  wts 

data(m(i),i=l,21)/28.9664d0,  28.964d0,  28.964d0,  28.964d0, 

1  28.964d0,  28.964d0,  28.962d0,  28.962d0, 

2  28.880d0,  28.560d0,  28.070d0,  26.920d0, 

3  26.660d0,  26.500d0,  25.850d0,  24.690d0, 

4  22.660d0,  19.940d0,  17.940d0,  16.840d0, 

5  16.170d0/ 
c  first  pass  flag 

data  ifirst  /  0  / 
c 

c  define  constants  on  first  pass 
c 

if(  ifirst  .ne.  0  )  go  to  1000 
c 

b=3.139d-7 

c  acceleration  of  gravity 
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g0=9.7803d0 

c  universal  gas  const,  j/kg 
rr=8313.432d0 
r=rr/m(l) 

c  planetary  radius,  meters 
re  =  6378145d0 
c 

c  initialize  lapse  rate  for  altitude  regions 
c 

do  10, 1=1,21 

lr(l)  =  ( tm(l+l)  -  tm(l) )/( z(l+l)  -  z(l) ) 

10  continue 
69  close(7) 

d0=p0/(r*tm(l)) 

p(l)=p0 

d(l)=d0 
do  20,1=1,20 
r=rr/m(l) 

if  (Ir(l).eq.O.dO)  then 

el=l  .d0-(b/2.d0)*(z(l+l)-z(l)) 
e2=g0*(z(l+l)-z(l))/(r*tm(l)) 
p(l+l)=p(l)*dexp(-el  *e2) 
d(l+l)=d(l)*dexp(-el  *e2) 
dd(l)=d(l+l)-d(l)/(z(l+l).z(l)) 

elseif  (lr(l).ne.0.d0)  then 

el=l.dO+(lr(l)/tm(l))*(z(l+l)-z(l)) 

e2=g0*b/(r*lr(l)) 

e3=e2*(z(l+l)-z(l)) 

e4=e2/b*(b/e2+l.d0+b*((tm(l)/lr(l))-z(l))) 
e5=e2/b*(LdO+  b*((tm(l)/lr(l))-z(l))) 
p(l+l)=p(l)*(el  **(“e5))*dexp(e3) 
d(l+l)=d(l)*(el  **(-e4))*dexp(e3) 
dd(l)=(d(l+l)-d(l))/(z(l+l).z(l)) 

endif 
20  continue 
c 

ifirst  =  1 
c 

1000  continue 
c 

c  determine  which  region  altitude  falls  into 
c 

do  500  j  =  1,20 

if(  alt  .It.  z(j+l) )  then 

go  to  501 
endif 

500  continue 
i  =  21 

501  continue 
c 

c  determine  parameters  at  altitude 
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talt=tm(i)+lr(i)*(alt-z(i)) 
galt=gO*(re*  *2/((re+alt)*  *2)) 
r=rr/m(i) 

if  lapse  rate  is  zero 

if  (Ir(i).eq.O.dO)  then 

el=LdO-(b/2.dO)*(alt-z(i)) 
e2=g0*  (alt-z(i))/ (r*  tm(i)) 
palt=p(i)*dexp(-l  .dO*el  *e2) 
dalt=d(i)*dexp(-  LdO*e  1  *e2) 
deldr  =  -b/2.d0 
de2dr  =  gO/(r*tm(i)) 

dddr  =  -d(i)*dexp(-el*e2)*(  el*de2dr  +  deldr*e2) 

if  lapse  rate  not  equal  to  zero 

elseif  (Ir(i).ne.O.dO)  then 

el=l.dO+(lr(i)/(tm(i)))*(alt-z(i)) 

e2=g0*b/(r*lr(i)) 

e3=e2*(alt-z(i)) 

e4=e2/b*(b/e2+l.dOH-b*((tm(i)/lr(i))-z(i))) 
e5=e2/b*(l  .d0+  b*((tm(i)/lr(i))-z(i))) 
palt=p(i)*(el  **(-e5))*dexp(e3) 
dalt=d(i)*(el  **(-e4))*dexp(e3) 
deldr  =  lr(i)/tm(i) 
de3dr  =  e2 

dddr  =  d(i)*(  -e4*(el**(-e4-LdO))*deldr 
1  +  (el **(-e4))*de3dr  )*dexp(  e3  ) 

endif 


return 


Appendix  B  Abort  Descent  Program 


This  section  contains  all  the  source  code  for  both  methods  of  descent.  The  Abortl  program 
models  the  constant  controlled  descent  while  Abort2  covers  the  heading  angle  controlled  descent.  The  only 
subroutines  not  included  are  Haming  and  Atm  which  are  also  used  in  the  Gravtum  program  in  Appendix  A. 
The  subroutine  Dynam  although  included  in  this  section  only  once  has  a  portion  of  the  code  that  was  only 
used  in  the  Abort2  program. 


Program  Abortl 

implicit  Real*8  (a-h,o-z) 
c 

c  3D  reentry  over  a  spherical  rotating  earth.  Vinhpp21-27 
c  assuming  controls  are  pitch  and  roll  which  are  held  constant  for  each  run 
c  but  cycled  through  all  possible  angular  configurations 
c 

c  input:  state  is  radius  r,  longitude,  latitude,  speed,  flight  path 
c  angle,  and  heading  angle 
c 

c  output:  radius,  longitude,  latitude,  speed,  time 
common/ham/t,x(230,4),f(230,4),err(230),n,h 
c 

common/param/Lalt,pitch,roll,pO,mu,R,M,S,Omega,MAC 
real*8  Lalt,t,h,x,f,err,pitch,roll,pO,mu,R,M,S 
real*8  Omega,  MAC 
integer  nxt,n 

c  Files  to  store  the  abort  landing  Zones 

open(unit=  13,file='c:\afit\landl  .out',status— unknown') 
open(unit=  14,file='c:\afit\land2.out',status='unknown') 
open(unit=  1 5,file='c:\afit\land3.out',status='unknown') 
open(unit=  1 6,file='c:\afit\land4.out',status='unknown') 
open(unit=  1 7,file='c:\afit\land5.out',status— unknown') 
open(unit=  1 8,file='c:\afit\land6.out',status='unknown') 

c  Files  to  store  the  velocity  data  at  nominal  landing 

open(unit=  1 9,file='c:\afit\Vel  1  .out',status=' unknown') 
open(unit=  20,file='c:\afit\Vel2.out',status='unknown') 
open(iuiit=  2 1  ,file='c:\afit\Vel3  .out',status=' unknown') 
open(unit=  22,file='c:\afit\Vel4.out',status='unknown') 
open(unit=23,file='c:\afit\Vel5.out',status-unknown') 
open(unit=24,file='c:\afit\Vel6.ouf,status-unknown') 
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c  constants:  (in  SI  Units) 
mu  =  3.98601dl4 
R=  6378.135d3 
p0-  101325.0d0 
M=  108864.0d0 
S  -  249,90d0 

Omega  =  0.000072921 158560d0 
MAC  =  12.0580d0 
Lalt  =  25000 

c  This  is  the  loop  which  cycles  through  all  angles  by  an  increment  of  5  degrees 

dol000i  =  30,  90,5 
pitch  =  i/57.29578 
write(*,*)  pitch 
do  2000j=  0,90,5 
roll-j/57.29578 
write(*,*)  roll 

c  open  each  state  vector  produced  by  gravity  turn  program  one  at  a  time 

open  (unit=4,  file  =  'c:\afit\abort Lout', status  =’unknown') 
c  open  (imit=8,  file  ==  'c:\afit\abort2.out', status  ^'unknown') 

c  open  (unit=9,  file  =  'c:\afit\abort3  .out', status  -  unknown') 

c  open  (unit=10,  file  =  'c:\afit\abort4.out', status  -unknown') 

c  open  (unit=l  1 ,  file  =  'c:\afit\abort5.out', status  -  unknown') 

c  open  (unit=12,  file  =  'c:\afit\abort6.out’,status  -  unknown') 

c  read  in  state  from  gravity  turn  program  to  be  used  as  initial  state 

c  (radius,longitude,  latitude, velocity, flight  path  angle,  heading  angle) 

read(4,*)x(l,l),x(2,l),x(3,l),x(4,l),x(5,l),x(6,l) 

x(2,l)-x(2,l)/57.29578 

x(3,l)=x(3,l)/57.29578 

n  =  6 

h  =  .05d0 

t  =  0.d0  . 

nxt  =  0 

call  haming(nxt) 
if(nxt.eq.0)then 
write(*,*)  'dead' 
stop 
endif 

do  while(x(l,nxt)-R.ge.Lalt) 
call  haming(nxt) 

c  write  to  the  correct  file  for  each  abort  run 

write(13,*)  x(2,nxt)*57.29578,x(3,nxt)*57.29578 
write(19,*)  x(4,nxt),t 
endif 
close(4) 
end  do 

2000  continue 
1000  continue 
stop 
end 
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Program  Abort2 

implicit  Real*8  (a-h,o-z) 
c 

c  3D  reentry  over  a  spherical  rotating  earth.  Vinhpp  21-27 
c  assuming  controls  are  45  degrees  of  pitch  and  roll  until  the  proper  heading 
c  angle  to  intercept  the  landing  site  is  achieved  then  pitch  =  30  and  roll  =  0  degrees 
c 

c  input;  state  is  radius  r,  longitude,  latitude,  speed,  flight  path 
c  angle,  and  heading  angle 
c 

c  output:  radius,  longitude,  latitude,  speed,  time 

conunon/ham/t,x(230,4),f(230,4),err(230),n,hc 

common/param/Lalt,pitch,roll,pO,mu,R,M,S,Omega,MAC, 

1  LatSite,  Lon_Site 
real*8  Lalt,t,h,x,f,err,pitch,roll,pO,mu,R,M,S 
real*8  Omega,  MAC,  Lat_Site,  Lon_Site 
real‘''8  Lon_Plus,Lon_Minus,Lat_Plus,Lat_Minus 
integer  nxt,n,count 

open(unit=  8,file='c:\afit\contrl.out',status-unknown') 
open(unit=  1  l,file='c:\afit\AltVel.out',status='unknovvn') 


c  constants:  (in  SI  Units) 
mu  =  3.98601dl4 
R=  6378.135d3 
pO  =  101325.0d0 
M=  108864.0d0 
S  =  249.90d0 

Omega  =  0.000072921 158560d0 
MAC  =  12.0580d0 

c  Minimum  Altitude  over  center  of  runway  for  safe  landing 

c  estimated  at  20,000  ft  or  6098  meters 

Lalt  =  6098 
coimt=0 , 

c  open  state  vector  produced  by  gravity  turn  program 

open(unit=l  0,file  -c:\afit\grav.out',status- unknown') 

do  while(count.le.l80) 
count  =  coimt  +  1 

c  read  in  state  from  gravity  turn  program  to  be  used  as  initial  state 

c  (radius,longitude,  latitude, velocity,flight  path  angle,  heading  angle) 

read(10,*)x(l,l),x(2,l),x(3,l),x(4,l),x(5,l), 

1  x(6,l) 

x(2,l)=x(2,l)/57.29578 
x(3,l)=  x(3,l)/57.29578 

c  Input  the  longitude  and  latitude  of  the  east  coast  airport  chosen  for  landing 

c 

Lon_Site  =  -60.0/57.29578 
Lat_Site  =  46.17/57.29578 
c 
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Zone  for  tracking  intercept  of  airport  longitude  and  latitude  by  +/-  .05 

Lon_Plus  =  Lon_Site  +  .05/57.29578 
Lon_Minus  =  Lon_Site  -  .05/57.29578 
Lat^Plus  =  Lat^Site  +  .05/57.29578 
Lat_Minus  =  Lat_Site  -  .05/57.29578 

for  45  degree  pitch  and  45  degree  roll 

pitch  =  0.78540d0 
roll  =  0.78540d0 
n  =  6 
h  =  .05d0 
t  =  0.d0 
nxt  =  0 

call  haming(nxt) 
if(nxt.eq.0)then 
write(*,*)  'dead’ 
stop 
endif 

do  while(x(l,nxt)-R.gt.Lalt) 
call  haining(nxt) 

if(dabs(t-INT(t+.01)).lt..025)then 
write(*,*)  count 

write(8,*)  x(2,nxt)*57.29578,x(3,nxt)*57.29578 
if(x(2,nxt).gt.Lon_Minus.and.x(2,nxt) 
.lt.LonPlus.or.x(3,nxt).gt.Lat_Minus.and. 
x(3,nxt).lt.Lat_Plus)then 

write(ll,*)  x(2,nxt)*57.29578,x(3,nxt)*57.29578,x(l,nxt>R 
,x(4,nxt),t,count 
if(x(  1  ,nxt).le.Lalt)then 
write(l  1  ,*)  x(l  ,nxt)-R,  x(4,nxt),t,count 
if(x(2,nxt)*557.29578.gt.-65)then 
stop 
endif 
endif 
endif 
endif 
end  do 
end  do 
close(lO) 
stop 


Subroutine  Dynam(nxt) 
c 

c  3D  reentry  over  a  spherical  rotating  earth.  Vinh  pp  21-27 
c  calculates  eom  assuming  controls 
c  are  angle  of  attack  and  roll 
c 

c  state  is  radius  r,  longitude,  latitude,  speed,  flight  path 
c  angle,  and  heading  angle 
c 

implicit  Real*8  (a-h,o-z) 

common  /dyn/g,icl,clnew,clmax,rho,cdmax, 

1  aoamax,icount 

common  /param/  Lalt,pitch,roll,pO,mu,R,M,S,Omega,MAC, 
1  Lat_Site,  Lon_Site 

common  /ham/  t,x(230,4),f(230,4),err(230),n,h,mode 

double  precision  dynpre 

real*8  t,x,f,err,h,g 

double  precision  cd,cl 

integer  n,nxt 

double  precision  clmax,clnew,cdmax,aoamax 
double  precision  Lon  Site,  Lat_Site,  1,  Psi 
double  precision  Lalt, pitch, roll, pO,mu,R,M,S,Omega,  MAC 
double  precision  rho,alt 

c  constants 

mu  is  the  Gravitational  constant  of  the  Earth 
c  R  is  Radius  of  the  Earth 

c  POis  pressure  at  sea  level 

c  M  is  mass  of  vehicle 

c  S  is  vehicle  reference  area 

c  Omega  is  rotation  rate  of  the  Earth 

c  inverse  square  gravity 

g  =  mu/(  x(l,nxt)*x(l,nxt) ) 

c 

c  determine  controls  at  this  state 

c  call  contrl(nxt,cd,cl,dcddp,dcldp,rho,drhodr) 
alt  =  x(l,nxt)-R 
c 

call  atm(alt,pO,PALT,TALT,rho,dDdr,d2Ddr, sonic, mfjp,dm^dr) 
Kn  =  mfp/MAC 

call  aero(pitch,Kn,cd,cl,cdp,clp) 

c  dynamic  pressure  =  rho  S  v^2  /  2 
c 

dynpre  =  rho  *  S  *  x(4,nxt)*x(4,nxt)  /  2.d0 


B-5 


c  equations  of  motion 
c 

C  f(l ,nxt)  -  D(R)/DT  -  radius 
C  f(2,nxt)  -  D(theta)/DT  -  longitude 
C  f(3,nxt)  -  D^hi)/DT  -  latitude 
C  f(4,nxt)  -  D(vel)/DT  -  velocity 
C  f(5,nxt)  -  D(gamma)/DT  -  flight  path  angle 
C  f(6,nxt)  -  D(psi)/DT  -  heading  angle 

c  This  portion  is  only  used  in  the  Abort2  Program 
c  This  is  to  control  the  Orbiter  to  fly  straight  when  it  is  at  the 

c  proper  heading  angle  to  intercept  the  selected  airport 
c 

1  =  dacos(dcos(1.5708-x(3,nxt))*dcos(1.5708-Lat_Site) 

1  +  dsin(L5708“X(3,nxt))*dsin(l  .5708-Lat_Site) 

2  *  dcos(x(2,nxt)-Lon_Site)) 

Psi  =  dasin((dsin(1.5708-x(3,nxt))*dsin(x(2,nxt)-Lon_Site)) 

1  /dsin(l)) 

if(x(6,nxt).ge.(Psi+l  .5708))then 
roll  =  0.0d0 
pitch  =  0.5236d0 
endif 

c  Vinh  eqns  2-28 

f(l,nxt)  =  x(4,nxt)*dsin(x(5jnxt)) 

f(2^t)  =  x(4^t)*dcos(x(5,nxt))*dcos(x(6,nxt))/(x(l,nxt) 

1  *dcos(x(3,nxt))) 

f(3,nxt)  =  x(4,nxt)*dcos(x(5,nxt))*dsin(x(6,nxt))/x(l,nxt) 
c 

c  Vinh  eqns  2-3 1 ;  all  coriolis  terms  retained 

f(45nxt)  =  -  CD*dynpre/M  -  g*dsin(x(5,nxt)) 

f(5,nxt)  =  CL*dynpre*dcos(roll)/(M*x(4,nxt)) 

1  -(G  -  (x(4,nxt)*x(4,nxt)/x(  1  ,nxt)) ) 

2  *dcos(x(55nxt))/x(4,nxt) 

3  +  2.dO*Omega*dcos(x(3,nxt))*dcos(x(6,nxt)) 

f(6,nxt)  =  -x(4,nxt)*dcos(x(5,nxt))* 

1  dcos(x(6,nxt))*dtan(x(3,nxt))/x(  1  ,nxt) 

2  +  CL*dynpre*dsin(roll)/(M*dcos(x(5,nxt))*x(4,nxt)) 

3  +  2.dO*Omega*(  dtan(x(5,nxt))*dcos(x(3,nxt))*dsin(x(6,nxt)) 

4  -  dsin(x(3,nxt))) 
c  if(mode  .eq,  0)  return 


return 

end 
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Subroutine  Aero(  alpha,  Kn,  Cd,  Cl,  Cdp,  Clp ) 


c 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 

c 


aerodynamic  model  for  the  shuttle,  Blanchard  et  al,  JSR  31,  550 
converted  to  angle  of  attack  alpha  in  radians;  outputs  Cd,  Cl  not 
Ca,  Cn;  also  returns  first  derivatives. 

double  precision  alpha,Kn,Cd,  Cl,  Cdp,  Clp 
double  precision  Cnc,  Cac,  Cnf,  Caf,  Cncp,  Cacp,  Cnfp,  Cafp 
double  precision  Cnbar,  Cabar,alpha2,alpha3 
double  precision  Ca,Cn,Cap,Cnp,Cdf,Clf 

alpha2  =  alpha  *  alpha 
alphas  =  alpha2  *  alpha 

hypersonic  continuum 

Cnc  =  -0.839782d0  +  3.0012d0  ♦  alpha  -  0.303891d0  *  alpha2 
Cac  =  -0.0086314d0  +  0.190247d0  *  alpha  -  0.220613d0  *  alpha2 
1  +0.110351d0*  alphas 

Cncp  =  3.0012d0  -  0.60778 ld0*alpha 

Cacp  =  0.190247d0  -  0.441227d0  *  alpha  +  0.331053d0  ♦  alpha2 

free  molecular  flow 

Cnf  =  0.00158739d0  +  0.526217d0  *  alpha  +  3.17184d0  *  alpha2 
1  -L34772d0*  alphas 

Caf  =  0.751 105d0  +  0.944601d0  *  alpha  +  L94409d0  *  alpha2 
1  -2.20399d0*  alphas 

Cn£p  =  0.52621 7d0  +  6.34367d0  ♦  alpha  -  4.043 17d0  *  alpha2 
Cafp  =  0.944601d0  +  3.88819d0  *  alpha  -  6.61198d0  *  alpha2 

bridging  coefficients 

if(  dlogl0(  Kn  )  .It.  1.3849d0 )  then 
Cnbar  =  dexp(  -0.29981d0  *  (  1.3849d0  -  dlogl0(  Kn )  )♦* 

1  1.7128d0) 

else 

Cnbar  =  l.dO 


if(  dlogl0(  Kn )  .It.  1 .2042d0 )  then 
Cabar  =  dexp(  -0.2262d0  *  (  1.2042d0  -  dlogl0(  Kn )  )** 
1  1.8410d0) 

else 

Cabar  =  l.dO 
endif 
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merged  normal  and  axial  coefficients 

Cn  =  Cnc  +  ( Cnf  -  Cnc  )*Cnbar 
Ca  =  Cac  +  ( Caf  -  Cac  )*Cabar 

Cnp  =  Cncp  +  ( Cn^  -  Cncp  )  *  Cnbar 
Cap  =  Cacp  +  ( Cafp  -  Cacp  )  *  Cabar 

convert  to  Cl,  Cd 

Cd  =  Ca  *  dcos(alpha)  +  Cn  *  dsin(alpha) 

Cl  =  -Ca  *  dsin(alpha)  +  Cn  *  dcos(alpha) 

Cdp  =  Cap  *  dcos(alpha)  +  Cnp  *  dsin(alpha) 
1  -Ca  *  dsin(alpha)  +  Cn  *  dcos(alpha) 
Clp  =  -Cap*  dsin(alpha)  +  Cnp  *  dcos(alpha) 
1  -Ca  *  dcos(alpha)  -  Cn  *  dsin(alpha) 


return 
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